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ABSTRACT 


Computationally  Fast  Algorithms  for 
ASHA  Spectral  Estimation 

by 

Koji  Ogino  and  James  A.  Cadzow 

(ABSTRACT) 

The  high  performance  method  for  obtaining  an  AHMA  model 
spectral  estimate  of  a  wide-sense  stationary  time  series  has  been 
found  to  provide  typically  superior  performance  when  compared  to 
such  comtemporary  approaches  as  the  Box-Jenkins  and  nunrtnnim 
entropy  methods.  In  this  report,  fast  recursive  algorithmic 
implementations  of  the  high  performance  method  are  developed.  They 
are  recursive  in  the  sense  that  as  a  new  element  of  the  time  series 
is  observed,  the  parameters  characterizing  an  ARMA  spectral  estimate 
are  algorithmically  updated.  The  number  of  multiplications  and 
additions  required  at  each  recursive  stage  are  of  the  order  p  with 
p  being  the  number  of  denominator  coefficients  of  the  ABMA  model. 
Methods  of  modification  of  the  data  are  applied  to  achieve  a 
significant  computational  improvement.  The  development  is  predicated 
on  utilization  of  various  projection  operators. 


Chapter  1 


INTRODUCTION 


The  mathematical  development  of  digital  signal  analysis  has 
been  an  area  of  primary  concern  since  the  digital  computers  develop¬ 
ment  over  two  decades  ago.  The  analysis  of  the  frequency  character¬ 
istic  of  a  signal  is  of  particular  interest  in  the  field  known  as 
"time  series  analysis."  Time  series  analysis  encompasses  such 
areas  as  statistics,  economics,  and  communications.  Most  of  the 
work  in  time  series  analysis  has  been  carried  out  by  statisticians. 
More  recently,  however,  many  advancements  in  the  analysis  of  time 
series  have  been  made  in  the  field  of  signal  processing  based  on 
power  spectral  estimation  concepts  and  time  domain  analysis. 

The  need  for  power  spectral  estimates  arises  in  a  variety  of 
contexts,  including  the  measurement  of  noise  spectra  for  the  design 
of  optimal  linear  filters,  the  detection  of  narrow-band  signals  in 
wide-band  noise,  and  the  estimation  of  parameters  of  a  linear  system 
by  using  a  noisy  excitation. 

Current  methods  of  spectral  estimation  can  be  broadly  classified 
into  two  categories.  One  is  the  classical  approach  which  includes 
the  periodgram  method,  autocorrelation  methods  and  its  variants 
(Bartlett,  1953;  Blackman  and  Tukey,  1958;  Grenander  and  Rosenblatt, 
1957;  Jenkins  and  Watts,  1968;  Koopmans,  1974).  The  second  is 
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modem  power  spectral  density  estimations  based  on  parameters 
modeling.  This  includes  the  maximum  entropy  method  (Burg,  1967),  one- 
step  linear  prediction  (Parzen,  1969),  and  spectral  estimation  using 
ARMA  model  (Tretter  and  Steiglitz,  ly67;  Gutowskl,  Robinson  and 
Treitel,  1978).  In  practical  signal  processing  applications,  classical 
approachs  have  been  incorporated  by  many  researchers  and  users. 

This  is  because  classical  methods  are  fairly  easy  to  implement 
and  can  be  computed  efficiently  by  using  the  fast  Fourier  transform 
(Cooley  and  Tukey,  1965).  However,  the  spectral  estimates  obtained 
by  classical  methods  can  provide  unsatisfactory  results  when  the  data 
length  is  short.  For  example,  variance  of  estimates  is  large  and 
the  resolution  capability  of  noise  embedded  sinusoids  is  poor  in  such 
cases.  To  overcome  these  difficulties,  the  modem  spectral  estimation 
methods  were  developed.  These  methods  provide  better  spectral 
performance  than  classical  methods.  For  example,  one  of  the  widely 
used  modem  spectral  methods  referred  to  as  the  Maximum  Entropy 
method  (Burg,  1967)  possesses  better  resolution  capability  than  the 
classical  periodgram  approaches  for  short  data  lengths.  The  Maximum 
Entropy  method  is  classified  as  an  autoregressive  (AR)  model.  The  AR 
model  is  also  known  as  an  all-pole  model  which  uses  only  a  denominator 
polynomial  of  a  rational  model.  In  recognition  of  this  constraint, 
a  more  general  form,  the  autoregressive  and  moving  average  (ARMA) 
model  which  has  numerator  polynomials  as  well  as  denominator  poly¬ 
nomial  has  been  proposed.  A  variety  of  procedures  has  been  developed 
for  generating  ARMA  models.  One  of  these  methods  is  the  so-called  'high 
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performance'  ARMA  method  which  was  recently  developed  by  Cadzow  (1979). 
The  'high  performance'  ARMA  method  has  provided  excellent  spectral 
estimation  performance  when  compared  with  the  Maximum  Entropy  and 
its  variants.  However,  its  computational  efficiency  is  relatively 
burdensome . 

Recently,  attention  has  been  directed  towards  developing  'fast' 
spectral  estimation  algorithms.  These  include  the  generalized  Levinson's 
algorithm.  As  an  example,  it  is  possible  to  use  this  approach  for 
estimating  the  autoregressive  coefficients  of  a  p-th  order  AR  model 
with  the  number  of  required  additions  and  multiplications  being  on 
the  order  of  p  (i.e.,  0(p  )).  Recently,  Morf  developed  the  doubling 
algorithm  which  reduced  the  required  computations  to  0(p  log  p)  by 
using  the  divide  and  conquer  approach  (Morf,  1980).  More  recently, 
recursive  methods  which  have  an  ability  to  compute  necessary  parameters 
at  the  arrival  of  each  new  data  point  has  been  proposed  (Lee  and 
Morf,  1980).  This  algorithm  does  not  require  any  matrix  formulation 
and  the  computational  requirements  can  be  reduced  to  0(p)  to  update 
the  AR  model  parameters  with  each  new  data  sample. 

In  this  report,  the  development  of  fast  algorithms  for  the  high 
performance  spectral  estimation  method  is  treated.  To  begin  our 
development,  in  Chapter  2,  the  mathematical  definition  of  power 
spectral  density  function  is  stated  and  two  classical  methods  referred 
to  as  the  periodgram  and  the  autocorrelation  method  are  discussed. 

The  common  weakness  of  these  classical  techniques  are  examined.  In 
Chapter  3,  a  standard  procedure  of  modern  spectral  estimation. 


namely,  the  rational  function  model  is  discussed.  Modem  spectral 
linear  estimators  can  be  classified  into  three  types  of  models: 

(i)  AR  (Autoregressive)  model,  (ii)  MA  (Moving  Average)  model,  and 
(iii)  ARMA  (Autoregressive  and  Moving  Average)  model.  It  is  widely 
known  that  the  ARMA  model  is  a  desired  form  from  a  parameter  parsimony 
viewpoint.  In  Chapter  4,  the  'high  performance'  ARMA  spectral  estima¬ 
tion  is  described.  Although  this  method  gives  excellent  spectral 
performance,  the  computational  requirements  are  relatively  burdensome. 
To  achieve  a  higher  degree  of  computational  efficiency,  fast  algorithms 
are  developed  in  Chapter  5  and  data  modification  methods  are  intro¬ 
duced.  In  Chapter  6,  a  recursive  algorithm  which  requires  0(p) 
computations  at  the  arrival  of  each  new  data  sample  is  developed. 
Development  of  this  algorithm  is  predicated  on  various  projection 
operator  decompositions. 


CONVENTIONAL  SPECTRAL  ESTIMATIONS 


2.1  Introduction 

The  spectral  density  function  is  mathematically  defined  in 
Section  2.2.  Conventional  spectral  estimation  techniques  have  been 
developed  based  on  the  Fourier  transform  relationship  between  the 
power  spectral  density  function  and  the  autocorrelation  sequence 
(Bartlett,  1953;  Blackman  and  Tukey,  1958;  Grenander  and  Rosenblatt, 
1957;  Jenkins  and  Watts,  1968;  Koopmanns,  1974).  For  example  Blackman 
and  Tukey  developed  an  autocorrelation  method  (Blackman  and  Tukey, 
1958)  which  includes  following  steps: 

(i)  Estimate  the  autocorrelation  sequence  from  the  observed 
data ; 

(ii)  Window  the  autocorrelation  estimate; 

(iii)  Fourier  transform  of  the  windowed  data  record. 

While  various  procedures  are  used  in  step  (i)  to  estimate  the  auto¬ 
correlation  function,  the  objective  is  usually  to  obtain  a  minimum 
bias  and  minimum  variance  estimate  of  the  true  autocorrelation 
sequence.  In  step  (ii) ,  windowing  is  used  to  reduce  the  bias  and  the 
variance  of  the  power  spectral  estimate.  However,  the  windowing 
process  decreases  the  resolution  of  the  power  spectral  estimate. 

This  autocorrelation  method  demonstrates  typical  weaknesses  of 
conventional  spectral  estimation  approaches.  Spectral  estimation 
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performance  had  not  been  improved  until  the  development  of  modern 
spectral  estimation  techniques. 

2.2  Definition  of  Power  Spectral  Density 

Let  us  consider  a  discrete  time  random  sequence  {x(n)}  with 
autocorrelation  sequence  (rx(m)}  defined  by 

rx(m)  »  E  [x(n  +  m)  x*(n)]  (2.2.1) 

where  E  and  *  denote  the  expected  value  and  complex  conjugate  operation, 
respectively.  We  will  denote  the  z-transform  of  {rx(m)}  by 

oo 

8,(2)  -  Z  r  (m)  z"m  (2.2.2) 

X  _  X 

m=a-oo 

The  associated  power  spectral  density  is  then  defined  to  be 

OO 

S  (u)  *  S(z) I  «  Z  r  (m)  e"jwm  (2.2.3) 

X  X  ]  UlJ  _  X 

2— e  msB-oo 

Applying  the  inverse  z-transform  to  eq.  (2.2.2),  we  have 

',<»>  -  SJ  f  sx(2)  2'”  T  (2-2-4) 

C 

where  C  is  a  simple  closed  contour  contained  within  the  region  of 
convergence  for  Sx(z).  If  C  is  chosen  to  be  the  unit  circle,  by 
making  the  change  of  variable  z=e^U,  we  derive  the  discrete  inverse 
Fourier  transform  relationship 
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rx(m)  =  ^  Sx(u>)  eja“  da> 


(2.2.5) 


The  variance  of  the  random  time  series  (x(n) }  is  equal  to  rx(0)  and 


can  be  expressed  by 


E{ |x(n) j  }  =  rx(0) 


2tt  _ 


Sx(w)  diu 


(2.2.6) 


It  follows  that  the  average  power  in  the  incremental  frequency  band 


u>0  _<  a)  <_  +  dw  (Tretter,  1976)  is  found  to  be 


P  (cj  )  »S  (iDn)  T. — 
xv  0'  xv  0^  2 it 


(2.2.7) 


As  shown  in  eq.  (2.2.6),  the  time  series  variance  is  equal  to  the 
total  power  of  the  signal  which  is  a  scalar  multiple  of  the  area 
under  the  curve  Sx(ui).  Observing  the  relation  between  expressions 
(2.2.6)  and  (2.2.7),  one  can  see  that  the  integral  over  the  incremental 
frequency  band  is  proportional  to  the  total  power  of  the  signal  in 
that  band.  For  these  reasons  the  function  Sx(u>)  is  called  the  power 
spectral  density. 

The  frequency  response  of  a  linear  shift-invariant  system  and 
the  frequency  domain  representation  of  a  discrete-time  signal  are 
essential  concepts  in  digital  signal  processing.  In  this  section  we 
describe  another  interpretation  of  the  power  spectral  density 
function  using  the  theory  of  linear  discrete-time  systems  for  the 
case  when  the  input  is  a  random  time  series  (Oppenheim  and  Schafer, 
1975).  Consider  a  stable  linear  shift-invariant  system  with  unit- 
sample  response  h(n) .  Let  e(n)  be  a  real  input  sequence  that  is  a 


B 


V 
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sample  sequence  of  a  wide-sense  stationary  discrete-time  random 
process.  Then  the  output  of  the  linear  system  is  a  sample  function 
of  a  random  process  related  to  the  input  process  by  the  linear 
transformation 


x(n)  *  E  h(n  -  k)  s(k) 
k*-“ 


(2.2.8) 


It  can  be  shown  that  if  the  input  is  stationary,  then  so  is  the  output. 
The  input  signal  may  be  partially  characterized  by  its  mean  and  its 
autocorrelation  function  r^ (m) ,  or  we  may  also  have  additional 
information  about  first  or  higher  order  probability  distributions. 

In  characterizing  the  output  random  process  (x(n) } ,  we  desire  similar 
information.  For  many  applications,  it  is  sufficient  to  characterize 
both  the  input  and  output  in  terms  of  simple  averages,  such  as  means, 
variances,  and  autocorrelations.  Therefore,  we  shall  derive  input- 
output  relationships  between  these  quantities.  Generally  we  consider 
zero  mean  processes  and  our  analysis  is  restricted  to  the  examination 
of  the  autocorrelation  sequence.  The  autocorrelation  function  of  the 
output  process  is  readily  shown  to  be  given  by 


r  (m)  -  E  r  (m  -  n)  E  h(k)  h(n  +  k)  (2.2.9) 

X  £  , 

Ila— 00 


To  characterize  the  response  of  a  linear  time-invariant  system  to  a 
discrete  time  input,  we  apply  the  z-transformation  to  expression 


(2.2.9)  to  yield 


Sx(z)  -  H(z)  H(z)  S_(z) 


(2.2.10) 
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where  H(z)  is  the  transfer  function  of  the  linear  shift-invariant 
system.  In  terms  of  the  power  spectral  density,  (2.2.10)  becomes 


Sx(w)  *  |H(e-’<J)|2  S£(u)) 


(2.2.11) 


where  the  impulse  response  (h(k)}  is  taken  to  be  a  real  sequence. 

2 

If  the  input  random  process  is  a  white  noise  with  variance  0£  ,  it 
follows  that 


S  (u>)  =  lH(e^“)  | 2  a  2 


(2.2.12) 


Relationship  (2.2.12)  is  extensively  used  in  analysis  concerned  with 
modern  spectral  estimation. 


2.3  Discrete  Fourier  Transform  Approach 


As  shown  in  Section  2.1,  the  power  spectral  density  and  auto¬ 
correlation  functions  are  related  by  the  discrete  Fourier  transform. 
Suppose  that  the  sequence  (x(n)}  is  a  wide-sense  stationary  random 
time  series  and  the  complete  knowledge  of  the  associated  autocorrela¬ 
tion  (rx(m)}  is  given,  the  spectral  density  can  be  simply  obtained  by 


S  (w)  =  I  r  (m)  e 
x  x 

m*-® 


(2.3.1) 


In  relevant  signal  processing  applications,  it  is  never  feasible 
to  measure  an  infinite  number  of  autocorrelation  sequence  elements 


(rx(m) } . 
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We  will  now  begin  to  examine  the  problem  of  estimating  Sx(oj) 
from  a  finite  observation  of  the  time  series  (x(n)}.  This  observation 
can  be  represented  by  a  set  of  N  contiguous  samples 

x(0) ,  x(l),  ...  ,  x(N-l)  (2.3.2) 

About  two  decades  ago,  spectral  estimates  had  been  mostly  accomplished 
by  the  periodgram  and  autocorrelation  methods. 

2.3.1  Periodgr am  Method 

To  include  an  additional  degree  of  flexibility,  suppose  that  the 
observed  sequence  is  modified  to  form  the  auxiliary  signal 

f(n)  =  w(n)  x(n)  0  <_  n  <_  N-l  (2.3.3) 

where  w(n)  *  0  for  n  <  0  and  n  M.  The  sequence  w(n)  is  frequently 
called  a  data  window.  The  sample  autocorrelation  function  for  the 
modified  observed  sequence  can  be  written  as 

00 

r,(n)  *  i  Z  f(k+n)  f(k)  (2.3.4a) 

f  N 

=  i  f (n)  *  f(-n)  (2.3.4b) 

N 

where  *  denotes  the  operation  of  convolution.  Denoting  the  z-transform 
of  r^(n)  and  f(n)  by  R^(z)  and  F(z),  respectively,  the  convolution 
and  time  reversal  theorems  yield  the  following  relationship 

Rf(z)  -  i  F(z)  F(z_1) 


(2.3.5) 
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Evaluating  this  expression  at  z  =  we  have 

Rf(ej“)  =  i  |  Sf(w)  !  2  (2.3.6) 

The  function  R^(e^w)  is  known  as  the  periodgram  of  (f(n)}.  Two 
decades  ago,  the  periodgram  method  became  popular  because  R^(e',u) 
could  be  computed  efficiently  by  using  the  fast  Fourier  transform 
(FFT  see  Cooley  and  Tukey,  1965). 


2.3.2  Autocorrelation  Method 

When  the  true  autocorrelation  function  r  (m)  is  unknown,  it  is 

x 

desired  to  calculate  an  estimate  of  the  autocorrelation  function. 
The  associated  spectral  estimate  can  then  be  obtained  by  taking  a 
Fourier  transform  of  this  autocorrelation  estimate  (Blackman  and 
Tuckey,  1959) .  Two  common  estimates 


and 


r  (m)  =  - — 
x  N-m 


N-m 

I  x(i)  x(i+m) 


N-l 


(2.3.7) 


*  1  N_m  * 
rx(m)  *  jj-  2  x(i)  x(i+m) 

1=1  m  =  0,  ...  ,  N-l 


(2.3.8) 


are  typically  used  for  estimating  the  autocorrelation  function. 
Applying  the  expected  value  operation  on  expression  (2.3.7),  we  obtain 


E  [r  (m)]  -  l  E  [x(i)  x(i+m)] 

x  N— m  ^  ^ 


(2.3.9a) 


*  r  (m) 
x 


(2.3.9b) 


The  autocorrelation  estimate  rx(m)  is  seen  to  be  an  unbiased  estimate. 
On  the  other  hand,  one  can  similarly  show  that  rx(m)  is  a  biased 
estimate.  Because  rx(m)  is  an  unbiased  estimate,  it  might  be  thought 

A  J 

r  (m)  is  the  better  estimate.  For  several  reasons,  however,  r  (m) 

X  X 

is  sometimes  preferable  to  r^(m) .  First,  the  biased  estimate  does 
not  violate  a  property  of  a  valid  autocorrelation  functions,  that  is 

r  (0)  >  | r(m) j  (2.3.10) 

while  the  unbiased  estimate  can  violate  this  property.  Second,  the 
biased  estimate  produces  a  nonnegative  spectral  estimate,  while  the 
unbiased  estimate  may  not  (Burg,  1975).  Third,  the  mean-square  error 
for  the  biased  estimate  is  less  than  that  for  the  unbiased  method 
(Jenkins  and  Watts,  1968).  And  finally,  Parzen  provides  an  argument 
in  favor  of  the  biased  estimate  by  claiming  that  rx(m)  has  less 
variance  than  rx(m)  (Parzen,  1974). 

Various  procedures  may  also  be  used  to  estimate  the  autocorrela¬ 
tion  function.  The  objective  of  these  procedures  is  usually  to  obtain 
a  minimum  variance  estimate  of  the  true  autocorrelation  function. 
Similarly,  the  estimate  of  the  autocorrelation  function  is  windowed 
to  reduce  the  bias  and  variance  of  the  power  spectral  estimate,  but 
increases  its  statistical  stability.  Various  window  functions  have 
been  used  which  are  generally  unrelated  to  the  data  or  the  random 
process  being  analyzed.  Both  the  finite  record  length  of  the  auto¬ 
correlation  function  estimate  and  the  windowing  process  applied  to 
the  autocorrelation  function  decreases  the  resolution  of  the  power 


spectral  estimate.  An  additional  disadvantage  of  windowing  is  that 
unless  one  performs  good  windowing,  excessive  side  lobes  may  be 
introduced  in  the  power  spectral  estimate.  Side  lobes  may  be  reduced 
by  employing  well  designed  windows  but  we  then  lose  spectral  resolution 
particularly  when  the  data  record  is  short. 

The  autocorrelation  method  and  its  variants  were  developed  to 
achieve  better  spectrum  estimate  performance  in  comparison  to  the 
periodgram  method.  As  indicated  above,  however,  the  autocorrelation 
method  has  still  several  disadvantages.  These  disadvantages  had  not 
been  overcome  until  the  development  of  modem  spectral  estimation 
techniques . 


Chapter  3 


MODERN  SPECTRAL  ESTIMATION 


3.1  Introduction 


One  of  the  most  widely  used  models  for  spectral  estimation  is 
the  rational  model.  The  stochastic  time  series  {x(n)}  is  said  to 
have  a  rational  power  spectrum  if  its  power  spectral  density  can  be 
expressed  in  the  form 


Sx(w)  ■  jH(e^)  |  2  a2 


(3.1.1) 


2 

where  a  is  a  positive  constant  and  the  characteristic  rational 
function 


joi.  b„  +  b.  e  + 
_  B(eJ  )  0  1 _ 


H(eJ“)  = 


A(ej“) 


+  b  e 

_ 3. 


-jqu 


1  +  a^  e""^w  + 


+  a  e 
P 


-jpw 


(3.1.2) 


is  composed  of  the  ratio  of  the  polynomials  A(e^UJ)  and  B(eJll))  which 
may  have  real  coefficients  and  the  zeros  of  A(eja>)  are  all  contained 
within  the  unit  circle.  The  rational  power  spectral  density  (3.1.1) 
is  said  to  have  order  (p,  q)  and  its  zeroes  and  poles  are  seen  to 
occur  in  reciprocal  complex  conjugate  pairs. 

A  particularly  convenient  interpretation  on  how  a  stochastic 
time  series  with  rational  spectrum  may  arise  follows  directly  from 
the  characteristic  rational  function.  This  entails  treating  the 
characteristic  rational  function  (3.1.2)  as  being  the  transfer  function 
of  a  causal,  time- invariant  linear  system.  It  then  follows  that  this 


I 


system  will  be  characterized  by  the  recursive  equation 

q  P 

x(n)  =  Z  b.  e(n-i)  -  I  a  x(n-i)  (3.1.3) 

i=0  1  i=l  1 

where  the  time  series  (e(n)}  and  (x(n) }  are  taken  to  be  the  excitation 

and  response  signals,  respectively.  It  has  been  shown  in  section  2.1 

that  when  the  excitation  time  series  {e(n)}  is  a  zero  mean  stationary 

2 

white  noise  time  series  with  variance  a  ,  then  the  power  spectral 
density  of  the  response  time  series  is  given  by  relationship  (3.1.1). 
Thus  a  stationary  random  time  series  with  rational  power  spectral 
density  can  be  interpreted  as  being  the  response  of  a  stable  causal, 
time-invariant  linear  system  to  a  white  noise  excitation. 

The  general  linear  system  (3.1.3)  is  commonly  referred  to  as  an 
autoregressive-moving  average  (AKMA)  model  in  the  spectral  estimation 

literature.  This  ARMA  model  is  said  to  be  of  order  (p,  q)  and  it 

gives  rise  to  the  rational  spectrum  (3.1.2)  which  possesses  zeroes 
as  well  as  poles.  The  AKMA  model  is  the  most  general  of  rational 
spectrum  models  possible  and  its  a^  and  b^  coefficients  uniquely 
characterize  the  spectrum. 

In  the  spectral  estimation  literature,  most  of  activity  has  been 
directed  towards  the  special  class  of  ARMA  models  known  as  auto¬ 
regressive  (AR)  models.  An  AR  model  is  one  in  which  the  numerator 
polynomial  B(e^u)  is  equal  to  the  constant  b^.  As  such,  the  AR 
model  is  also  referred  to  as  an  all-pole  model  since  its  transfer 


function  is  specified  by 
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H(ejw)  = 


A(eJ“) 


(3.1.4) 


This  all-pole  model  is  the  one  most  often  used  in  spectral  estimation. 

Another  subclass  of  rational  spectrum  models  which  has  received 
attention  is  the  so-called  moving  average  (MA)  model  as  characterized 
by  ACe^)  =  1.  The  transfer  function  of  a  MA  model  is  given  by 
BCe^40)  and  it  is  therefore  also  referred  to  as  an  all-zero  model. 

In  summary.  Table  3.1  shows  the  rational  spectrum  associated  with 
each  of  these  models. 


3.2  Moving  Average  Model 

Many  conventional  methods  of  spectral  estimation  are  classified 
as  MA  models.  For  example,  the  periodgram  and  correlation  methods 
which  have  been  discussed  in  Section  2.3  can  be.  described  in  terms 
of  a  MA  model.  Generally,  little  attention  has  been  focused  on  MA 
models.  Welch  has  introduced  (Welch,  1967),  however,  a  MA  model 
technique  which  is  particularly  applicable  to  the  direct  computation 
of  a  power  spectrum  estimate  that  uses  the  FFT.  In  this  technique, 
the  data  record  is  first  sectioned  into  K  =  N/M  segments  of  M  samples 
each  as  defined  by 

x(i)  (n)  =  x(n  +  iM  -  M)  0  <_  n  <_  M-l,  1  <  i  <  R  (3.2.1) 


A  window  w(n)  is  next  applied  directly  to  the  data  segments  before 
computation  of  the  periodgram.  Then,  the  K  modified  periodgrams  as 
specified  by 


Model 


Table  3. 1  Rational  Spectrum  Models 
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4i}  (a.)  -  ^  I  S  X(l)(a)  v(n)  e"jun  |  2  i  =  1,  2, 
n=0 


.  ,  K 
(3.2.2) 


are  computed,  where 


1  2 
U  *  —  Z  w  (n) 

M  n=0 


(3.2.3) 


and  the  final  spectrum  estimate  is  defined  as 


1  K  (i) 

B„(«)  -  -  Z  (<o) 

x  K  .„i  M 


(3.2.4) 


By  taking  average  of  periodgrams  of  each  data  segment,  the 
desired  smoothed  periodgram  is  obtained.  In  using  this  segmentation, 
the  variance  of  the  spectrum  is  reduced.  The  price  paid  for  this 
reduction,  however,  is  a  loss  in  frequency  resolution  and  an  increased 
bias  of  the  estimate. 


3.3  Autoregressive  (AR)  Model 

In  the  last  decade,  much  attention  has  been  focused  on  the 
analysis  of  AR  models.  Two  major  spectrum  estimation  methods  for  AR 
models,  referred  as  one-step  linear  prediction  and  the  maximum 
entropy  method  (MEM)  appeared  in  the  literature  of  mathematical 
statistics  (Parzen,  1969)  and  geosciences  (Burg,  1967;  Lacoss,  1971; 
Ulrvch,  1972).  Although  these  two  methods  take  different  approaches, 
it  has  been  shown  that  they  give  the  same  spectral  estimate  (A  van  den 
Bos,  1971). 
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3.3.1  One-Step  Linear  Prediction 

In  the  application  of  one-step  linear  prediction,  one  seeks  to 
characterize  the  spectral  density  of  a  time  series  based  upon  a  finite 
set  of  time  observation 

x(l),  x(2) ,  ...  ,  x(N)  (3. 3. 1.1) 

As  described  in  Section  3.1,  the  AR  model  is  structured  by 

x(n)  +  a^  x(n  -  1)  +  ...  +  a^  x(n  -  p)  =  e(n)  (3. 3. 1.2) 

in  which  e(n)  is  a  white  noise  time  series  with  zero  mean  and  variance 
2 

.  The  objective  of  spectral  estimation  will  be  that  of  modeling 
an  underlying  time  series  {x(n) }  with  the  AR  model  structure  (3. 3. 1.2) 
in  which  the  a^  coefficients  are  estimated  from  the  given  finite  set 
of  observations  (3. 3. 1.1).  This  is  readily  achieved  by  applying 
the  well  known  method  of  one-step  linear  prediction. 

A  p-th  order  one-step  linear  prediction,  by  definition,  estimates 
the  value  of  a  random  time  series  using  a  linear  combination  of 
the  most  recent  p  samples.  Namely,  the  sample  x(n)  is  estimated  by 
means  of  the  relationship 
P 

x(n)  =  -  E  a  x(n  -  k)  (3. 3. 1.3) 

k=l  k 

The  difference  between  this  predicted  value  and  the  observed  value 
x(n)  over  the  observation  interval  is  called  the  prediction  error 
and  is  specified  by 

e(n)  =  x(n)  -  x(n)  p  <  n  <_  N  (3. 3. 1.4) 
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or 


P 

e(n)  *  x(n)  +  Z  x(n  -  k)  p  <  n  <_  N 
k«l 

Writing  these  error  expressions  in  matrix  form  yields 


e  »  x  +  X  a 


(3. 3. 1.5) 


(3. 3. 1.6) 


whete  a,  e^  and  x  are  p  x  1,  (N  -  p)  x  1,  and  (N  -  p)  x  1  column 
vectors,  respectively,  given  by 


a  -  [a^  ...  ,  ap]T 
e  =  [e(p  +  1),  e(p  +  2) ,  . . .  ,  e(N)]X 


(3.3.1.7a) 

(3.3.1.7b) 


X  =  [x(p  +  1),  x(p  +  2),  ...  ,  x(N)]T 


(3.3.1.7c) 


and  X  is  an  (N  -  p)  x  p  matrix  specified  by 


x(p)  x(p  +  1) 
x(p  -  1)  x(p) 


x(N  -  1) 
x(N  -  2) 


T 


x ( 1)  x(2) 


x(N  -  p) 


(3.3.1.7d) 


v/here  the  superscript  T  denotes  the  transpose  operation. 

The  a^  coefficients  are  to  be  now  selected  so  as  to  cause  each 
of  the  prediction  error  terms  e(n)  to  be  close  to  zero.  This 
selection  process  will  give  rise  to  the  so-called  optimal  one-step 
predictor.  To  achieve  the  required  objective  of  setting  the  e(n)  to 
be  near  zero,  one  typically  appeals  to  the  least  squares  method  which 
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minimizes  a  squared  error  criterion  of  the  form 


f  (a)  =  e  We 


(3. 3.1.8) 


where  W  is  an  (N  -  p)  x  (N  -  p)  nonnegative  definite  square  matrix. 
The  minimization  of  this  quadratic  functional  with  respect  to  the 
column  vector  a  is  straightforwardly  carried  out  and  results  in 


XT  W  X  a0  »  XT  W  x 


(3. 3. 1.9) 


It  can  be  shown  that  the  resulting  power  spectral  density  estimate 
of  the  time  series  {x(n) }  is  then  given  by 


S  (u) 
x 


1  +  a?  e-jtU  +  a°  e‘j2“  +  ...  ^‘’e^2 

12  p  1 


(3.3.1.10) 


where  the  a^  coefficients  are  obtained  upon  solving  relationship 
(3. 3. 1.9). 


3.3.2  Maximum  Entropy  Method  (MEM) 

The  MEM  is  a  result  of  Burg's  attempt  (Burg,  1967)  to  derive 
a  procedure  for  increasing  spectral  resolution  when  only  a  small 
number  of  samples  or  estimates  of  autocorrelation  function  are  avail¬ 
able.  As  mentioned  in  Section  2.3.2,  in  the  autocorrelation  method 
one  first  estimates  the  autocorrelation  function,  append  zeroes  to  in¬ 
crease  the  length  of  the  estimated  autocorrelation,  and  then  applies  the 
Fourier  transform.  In  contrast,  the  MEM  suggests  that  the  estimated 
autocorrelation  function  should  be  extrapolated  beyond  the  data 
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limited  range.  The  principle  used  for  this  extrapolation  process  is 
that  the  spectral  estimate  must  be  the  most  random  or  have  the  maximum 
entropy  of  any  power  spectrum  which  is  consistent  with  the  sample 
values  of  the  estimated  autocorrelation. 

In  the  analysis  of  MEM,  it  is  assumed  that  we  possess  a  partial 
autocorrelation  sequence  (r(0),  r(+l),  ...  ,  r(+M)}  which  is  a  subset 
of  a  infinite  extent  autocorrelation  function  (r(0),  r(+l) ,  ......}. 

It  is  desired  that  we  produce  from  this  partial  autocorrelation 
sequence  a  spectral  representation 

CO 

Sr(u)  *  £  r(k)  e"j“n  (3. 3.2.1) 

IT*-” 

which  is  a  Fourier  transform  of  the  autocorrelation  function  of 
infinite  length.  For  some  spectral  density  function  S^(w),  we  may 
associate  a  time  series  (f(n)}  by  means  of  inverse  Fourier  transform 
ir 

S^(ui)  e^un  doj  for  n  =  0,  +1,  ...  (3. 3. 2. 2) 

-77 

so  that 


r(n)  =  f(n)  for  n  =  0»  +1,  ...  ,  +M  (3. 3. 2. 3) 

This  expression  does  not  provide  us  with  a  unique  expression  for  the 
3pectrum  Sr(oi) .  To  overcome  this  difficulty.  Burg  developed  a  new 
spectral  estimator  called  the  maximum  entropy  method  (Burg,  1967). 

The  entropy  associated  with  power  spectrum  density  Sr(u>)  is  defined 


to  be 
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H 


log  [Sr  <a>)  ]  da) 


(3. 3. 2. 4) 


-IT 


Maximizing  the  entropy  with  respect  to  the  unknown  r(n)  for  |n|  >  M 
with  the  constraint 


r(n) 


■&!  Sr 


(w)  e^11  dw  for  Ini  >  M 


(3. 3. 2. 5) 


-TT 


results  in  the  maximum  entropy  spectral  estimate.  This  estimate 
expresses  maximum  uncertainty  with  respect  to  the  unknown  information 
that  is  consistent  with  the  known  information.  The  problem  of 
estimating  Sr(w)  becomes  a  calculus  of  variations  problem.  The  solution 
procedure  which  begins  with  the  introduction  of  a  Lagrange  multiplier 
for  each  of  the  constraint  equations  is  not  difficult  and  results 
in  the  spectral  estimate  (Burg,  1967) 


Sr(u)  = 


|l  +  a°e^U  +  —  + 


o  -JMco.2 


(3. 3.2.6) 


where  optimum  selection  of  coefficients  (k  *  1,  ...  ,  M)  are 
obtained  by  solving  the  following  matrix  system  of  equations 


r(0) 

r(l) 

.  .  .  r(M) 

r(l) 

r(0) 

.  .  .  r(M-l) 

r(M) 

r(M-l) 

.  .  .  r(0) 

(3. 3.2. 7) 


24 


Equation  (3. 3.2. 7)  can  be  solved  efficiently  by  using  Levinson’s 

2 

Algorithm  which  requires  0(M  )  computations  (Levinson,  1947). 

3.4  ARMA  Model 

A  variety  of  procedures  have  been  developed  for  generating  AHMA 
spectral  models.  These  include  the  whitening  filter  approach  which 
is  typically  iterative  in  nature,  generally  slow  in  convergence,  and, 
usually  requires  an  excessively  large  number  of  time  series'  obser¬ 
vations  to  be  effective  (Tretter  and  Steiglitz,  1967;  Gutowski, 

Robinson  and  Treitel,  1978) .  More  desirable  closed  form  procedures 
which  overcome  these  deficiencies  have  been  offered.  These  include 
the  so-called  Box-Jenkins  method  and  its  variants  (Box  and  Jenkins, 
1976;  Kaveh,  1979;  Kinkel,  Perl,  Scharf  and  Stubberud,  1979),  and, 
more  recently,  Cadzow  has  developed  a  "high  performance"  method 
(Cadzow,  1981).  In  this  section,  three  ARMA  methods,  namely,  the 
Whitening  method,  Gutowski  ARMA  method  and  Box-Jenkins  method  are 
briefly  discussed. 

3.4.1  Whitening  Method 

If  we  assume  that  the  Gaussian  random  series  (x(n)}  is  given, 
the  method  of  maximum  likelihood  (Haykin,  1979)  can  be  used  to  estimate 
the  coefficients  of  rational  spectrum  in  the  following  way.  Suppose 
the  time  sequence  (x(n)}  is  passed  through  a  transfer  function 
A(e^u)/B(e^a))  to  give  the  output  sequence  (e(n)}.  The  spectrum  of 


{ e (n) }  is  given  by 


If  one  could  choose  the  coefficients  of  A(e^W)  and  B(e^W)  so  that 
2 

S^Cu)  =  a ,  the  spectral  density  of  (x(n)}  would  be  given  by 


S  (u>) 
x 


A(ej“) 


(3. 4. 1.2) 


In  this  case,  (e(n)}  is  a  white  Gaussian  process.  The  maximum  likeli¬ 
hood  parameter  estimation  is  equivalent  to  finding  the  minimum  of  a 
function  of  several  variables  (Tretter  and  Steiglitz,  1967).  This 
is  called  the  minimum  residual  criterion  and,  intuitively,  one  attempts 
to  "whiten"  {x(n)}  as  much  as  possible.  The  whitening  process  is  sugges¬ 
tively  depicted  in  Fig.  3.4. 1.1. 

Because  of  the  rational  spectrum  model's  structure,  the  minimum 
residual  criterion  leads  to  nonlinear  equations  which  cannot  be 
solved  explicitly.  This  suggests  the  using  of  an  iterative  technique 
to  optimize  the  denominator  and  numerator  coefficients.  Many  such 
techniques  are  available,  ranging  from  steepest  descent  to  the 
Newton-Raphson  algorithm. 


3.4.2  Gutowski  ARMA  Method 

This  section  discusses  the  theoretical  motivation  for  the  ARMA 
modeling  technique  described  by  Gutowski  (Gutowski,  Robinson, 
Treitel,  1978).  Consider  the  discrete  time  linear  system 


shown  in  Fig.  3. 4. 2.1  with  input  u(k) ,  output  x(k) ,  and 
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Impulse  response  h(k).  If  the  transfer  function  H(z)  is  assumed  to  be 
a  rational  function  of  z,  then  it  may  be  written  as 

H(z)  -  (3. 4. 2.1) 

A(z) 

where  A(z)  and  B(z)  are  polynomials  of  z of  order  p  and  q,  respectively. 
This  assumption  in  turn  implies  that  the  output  is  described  by 

X(z)  =  u(z)  (3. 4. 2. 2) 

A(z) 

where  X(z)  and  U(z)  denote  z-transform  of  {x(k) }  and  (u(k)),  respectively. 
Gutowski's  ABMA  method  assumes  that  u(k)  is  equal  to  the  Kronecker 
delta  function  and  it  therefore  follows  that 


SK--xw  (3-4-2-3> 

Gutowski's  method  uses  Equation  (3. 4. 2. 3)  in  an  iterative  procedure 
to  estimate  A(z)  and  B(z)  from  the  data  sequence  (x(k) }.  Each  iteration 
may  be  described  in  terms  of  the  following  three  equations: 


(3. 4. 2. 4) 

(3. 4. 2. 5) 

(3. 4. 2. 6) 


The  basic  iterative  technique  may  be  seen  by  using  equation  (3. 4. 2. 4) 
through  (3. 4. 2. 6)  and  assuming  that  one  starts  with  a  reasonably  good 
estimate  of  B(z).  At  k-th  iteration,  the  following  steps  are 
required. 


(k)  (k)  (k) 

Step  1  Compute  A(z)  with  X(z)  input  and  B(z)  as  desired  output. 

(k)  (k) 

Step  2  Compute  C(z)  by  synthetic  division  of  the  value  1  by  A(k) . 

(k)  (k)  (k) 

Step  3  Compute  B(z)  with  C(z)  as  input  and  X(z)  as  desired 

output. 

(k)  (k) 

After  each  iteration,  if  A(z)  and  B(z)  are  better  than  the  previous 
iteration,  then  the  fit  will  improve.  At  the  completion  of  m-th 
iterations,  the  ARMA  spectral  estimate  is  given  by 


Sx(co) 


,  B(m)(ej%  2 
A(m)(ej“) 


(3. 4.  2. 7) 


The  above  procedure  is  repeated  until  convergence  occurs.  The 
minimum  delay  characteristics  of  Am(z)  is  guaranteed  by  the  fact  that 
the  inverse  is  computed  using  a  Toeplitz  formulation.  This  is  the 
strong  point  of  this  algorithm. 


3.4,3  Box-Jenkins  Method 

The  ARMA  model  with  order  (p,  q)  can  be  characterized  by  the 
following  recursive  relationship 

p  q 

x(n)  *  -E  a.  x(n-k)  +  l  b.  e(n-k)  (3. 4. 3.1) 

k-1  k=0  k 

n  *  p  +  1,  ...  ,  +  ® 

2 

where  (e(k))  is  a  white  noise  with  variance  a  .  The  autocorrelation 

G 


function  of  the  mixed  process  may  be  derived  by  multiplying  each 


(3. 4. 3. 2) 


side  of  (3. 4. 3.1)  by  x  (n-m)  and  taking  expectations  to  yield 


p  q 

rx  *  1  \  rx  (m  "  k)  +  E  bk  rX£  (m  "  k) 

k=l  k=0 


where  (n)  and  rx£(n)  denote  the  autocorrelation  of  the  sequence 
{x(k)}  and  cross  covariance  function  between  (x(k)}  and  {e(k)}, 
respectively.  Since  x(n-k)  depends  only  on  inputs  which  have  occurred 
up  to  time  n-k,  it  then  follows  that 


rx£(n)  =0  n  >  0 


rx£(n)  *  0  n  <  0 


We  see  that  (3. 4. 3. 2)  implies 


P 

r  (n)  =  -  E  a,  r  (n-k)  for  n  ^  q  +  1 

X  ,  -  tC  X 

k=l 


and  yields  the  following  matrix  system  of  equations 


-1 

" 

- 

mm.  — 

(q)  ...  rx  (q-p+1) 

al 

rx  (t*+1) 

•  •  • 

•  mm 

rx  (q+p-1)  ...  rx  (q) 

a 

P 

It 

rx  (q+p) 

. 

(3.4.3.3a) 

(3.4.3.3b) 


(3. 4. 3. 4) 


(3. 4. 3. 5) 


The  coefficients  will  be  obtained  by  solving  the  equation  (3. 4. 3. 5). 
The  numerator  dynamics  of  the  ARMA  model  is  characterized  by  c^ 
coefficients  (Kaveh,  1979)  which  can  be  expressed  as 


Chapter  4 


HIGH  PERFORMANCE  ARMA  MODEL 

4.1  Introduction 

It  is  widely  recognized  that  an  ARMA  spectral  model  is  generally 
the  most  effective  linear  rational  model  from  a  parameter  parsimony 
viewpoint  (see  Section  3.1).  In  recognition  of  this  fact,  a  variety 
of  procedures  have  been  developed  for  generating  ARMA  models 
(Steiglitz,  1977;  Box  and  Jenkins,  1976;  Kaveh,  1979;  Kinkel,  Perl, 
Scharf  and  Stubberud,  1979).  Some  of  these  methods  were  discussed  in 
Section  3.4.  As  indicated  in  Section  3.4,  it  is  recognized  that 
these  methods  share  certain  deficiencies.  To  overcome  these 
deficiencies,  the  'high  performance'  ARMA  method  was  developed 
(Cadzow,  1979,  1980,  a,b) .  It  provides  an  excellent  spectral  estima¬ 
tion  performance  when  compared  with  other  spectral  estimation  methods. 
In  this  chapter,  the  'high  performance'  method  is  described  and 
numbers  of  numerical  examples  are  provided.  This  chapter  is  basically 
identical  to  references  (Cadzow,  1979,  1980  a,  b) .  The  development  of 
this  method  is  based  upon  some  fundamental  concept  governing  ARMA 
time  series  which  will  be  discussed  in  next  section. 


4.2  Fundamental  Concepts 


The  stationary  random  time  series  (x^l  whose  power  spectrum  is 
of  a  rational  form  may  be  modeled  as  the  response  of  the  causal  ARMA 
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system  of  order  (p,  q) 


x 


k-i 


e 


k-i 


(4.2.1) 


where  the  time  series  { e.  }  is  taken  to  be  a  zero  mean  white  noise 

k 

excitation  signal.  The  autocorrelation  description  of  this  system 

is  obtained  by  first  multiplying  each  side  of  expression  (4.2.1)  by 
* 

the  entity  x^  m  and  then  taking  the  expected  value.  This  results  in 
the  well  known  Yule-Walker  equations  as  specified  by 
P 

r  (m)  +  E  a.  r  (m  -  i)  =  0  for  m  >  q  (4.2.2) 

x  ,  _  i  x 

1=1 


The  Yule-Walker  equations  (4.2.2)  will  serve  as  the  basis  for  esti¬ 
mating  the  ARMA  model's  denominator  coefficients  (i.e.,  coefficients). 


4.3  Denominator  Coefficient  Selection 

In  this  section,  a  novel  procedure  for  estimating  an  ARMA  model's 
denominator  coefficients  shall  be  presented  (Cadzow,  1979,  1980  a). 
This  development  is  begun  by  first  evaluating  the  model  equation 
(4.2.1)  over  the  integer  set  p  +  1  <_  k  <_  n  to  obtain  the  time  series 
relationships 
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p+l  p  P-q+1 

*  bo' 

p+2  p+1  p-q+2 

bl 

£  £  ,  .  .  .  £ 

b 

n  n-1  n-q 

q 

It  will  be  compactly  written  in  the  matrix  format 
x  +  X  a_  =>  £  _b 


(4.3.1) 


(4.3.2) 


where  x,  a^  and  j)  is  (n-p)xl,  pxl  and  (n-p)xl  column  vector,  respectively. 
The  symbols  X  and  &  denote  (n-p)xp  and  (n-p)  x  (q+1)  Toeplitz  type 
matrices,  respectively.  The  entries  of  these  vectors  and  matrices 
are  directly  obtained  from  expression  (4.3.1). 

It  is  now  desired  to  utilize  relationship  (4.3.1)  in  conjunction 
with  the  Yule-Walker  equations  (4.2.2)  to  effect  a  procedure  for 
estimating  the  ARMA  model's  autoregressive  coefficients.  As  we  will 
see,  this  objective  is  attained  by  first  introducing  the  following 
(n-p)xt  Toeplitz  type  matrix 


p-q 


p-q-1 


V  ,  ,  A 

p-q+1  p-q 


x  ,  x  _ 
n-q-1  n-o-2 


p-q-t+1 


p-q-t+2 


n-q-t 


(4.3.3) 
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where  the  convention  is  adopted  of  setting  to  zero  any  matrix  entry 
x^  for  which  k  lies  outside  the  observation  set  1  £  k  £  n.  The 
integer  t  which  specifies  the  number  of  columns  of  this  matrix  will 
also  be  found  to  correspond  to  the  number  of  Yule-Walker  equations 
that  are  being  approximated  (i.e.,  relationship  (4.2.2)  for 
q  <  m  _<  q  +  t) .  It  thus  follows  that  this  integer  parameter  must 
be  selected  to  at  least  equal  p  (i.e.,  t  ^  p)  so  as  to  assure  a 
well  defined  set  of  equations  for  the  p  autoregressive  coefficients. 

The  above  mentioned  Yule-Walker  equation  approximation  is 
achieved  by  premultiplying  each  side  of  relationship  (4.3.2)  by  the 
complex  conjugate  transpose  of  matrix  Y  as  denoted  by  Y  to  yield 

Yr  x  +  YX  X  a  =  T  €  a  (4.3.4) 


To  demonstrate  that  this  system  of  equations  yields  a  logical  choice 
for  the  Yule-Walker  equation  approximations,  let  us  now  take  the 
expected  value  of  each  of  its  sides.  This  is  found  to  result  in 


P 

(n  -  m)  (r  (m)  +  Z  a.  r  (m  -  k) }  =■  0  (4.3.5) 

k=l  k  X 

for  q  <  m  <_  q  +  t 


Thus,  the  system  of  linear  equations  (4.3.4)  is  seen  to  provide  an 

unbiased  estimate  of  the  underlying  Yule-Walker  equations.  It  is  to 

be  noted  that  the  right  hand  side  term  has  zero  expected  value  due  to 

the  fact  that  the  expected  value  of  the  matrix  Y  £  is  the  null 

matrix.  This  is  a  direct  consequence  of  the  ARMA  model's  causality  and 

the  whiteness  of  the  excitation  process  which  results  in  E  {x  e,}  *  0 

n  k 
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for  all  n  <  k. 

With  these  thoughts  in  mind,  a  logical  procedure  for  selecting  the 
ARMA  model's  autoregressive  coefficients  is  suggested.  Namely,  they 
will  be  selected  so  as  to  cause  the  left  hand  side  of  relationship 
(4.3.4)  to  be  close  to  its  expected  value  which  is  the  zero  vector 

X 

(i.e.,  E  {Y  b }  =  0) .  If  this  selection  procedure  is  adopted,  an 
approximation  of  the  Yule-Walker  equations  which  in  some  sense  is 
"most  consistent"  with  the  given  time  series  observations  is  at  hand. 

A  computationally  tractable  measure  of  the  closeness  to  which  the 
left  side  of  relationship  (4.3.4)  is  to  the  zero  vector  is  provided 
by  the  following  quadratic  functional 

f (a)  =  [Y+  x  +  Yf  X  a]TA  [Y~  x  +  Y^  X  a]  (4.3.6) 

in  which  A  is  a  t  x  t  positive-semidef inite  diagonal  matrix  whose 
diagonal  elements  are  chosen  to  possibly  weight  differently  various 

X  X 

elements  of  the  error  vector  Y  x  +  Y  X  a.  It  is  a  simple  matter  to 
show  that  a  minimizing  autoregressive  coefficient  vector  must  satisfy 
the  consistent  system  of  p  linear  equations 

Xf  Y  A  Y+  X  a  -  -X1"  Y  A  Y7  x  (4.3.7) 

in  the  p  autoregressive  coefficient  unknowns.  One  then  solves  this 
system  of  p  equations  for  the  most  data  consistent  set  of  auto¬ 
regressive  coefficient  estimates. 
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4.4  Numerator  Dynamics 

A  variety  of  procedures  exists  for  determining  the  numerator 
dynamics  of  an  ARMA  time  series  once  the  AR  coefficients  have  been 
estimated.  In  this  section,  two  procedures  which  have  been  found  to 
be  particularly  effective  shall  be  described.  Each  makes  use  of  th.^ 
governing  ARMA  relationship  that  models  the  underlying  time  series. 

4.4.1  Yule-Walker  Equation  Method  (Cadzow,  1979) 

In  this  approach  to  estimating  the  numerator  dynamics,  we  first 
introduce  the  so-called  causal  image  of  a  time  series  autocorrelation 
sequence  as  specified  by 

r*(n)  =  -h  rx(0)  <5(n)  +  rx(n)  u(n)  (4.4. 1.1) 

in  which  S(n)  and  u(n)  designate  the  unit-sample  and  unit-step 
sequences,  respectively.  Making  use  of  the  complex  conjugate 
symmetrical  property  of  stationary  autocorrelation  sequences,  it  then 
follows  that  the  autocorrelation  sequence  can  be  uniquely  recovered 
from  its  causal  image  according  to  the  simple  relationship 

rx(n)  *  r+(n)  +  rx  (-tO*  (4. 4. 1.2) 

Upon  taking  the  discrete-Fourier  transform  of  this  relationship, 
it  follows  that  the  time  series  spectral  density  is  given  by 

Sx(a>)  -  S+(«)  +  [Sx(w)]* 

*  2  R  [S+(ui)]  (4. 4. 1.3) 

6  X 


where  S  (to)  denotes  the  discrete  Fourier  transform  of  the  causal 
x 

image  r+(n) .  According  to  relationship  (4. 4. 1.3),  one  may  attain  a 
spectral  density  estimate  by  estimating  S+Cw) .  This  will  be  the 
approach  taken  in  this  section. 

An  estimation  of  the  Yule-Walker  equations  (4.2.2)  which  govern 
the  ARMA  model  time  series  indicates  that  the  causal  image  sequence 
will  generate  the  auxiliary  (c^)  sequence  according  to 


c 

m 


=  r+<»>  + 


P 

E 

k-1 


\ 


+ 

r 

x 


(m  -  k) 


(4. 4. 1.4) 


m  =  0,  1,  ...  ,  s  for  s  =  max  (q,p) 


It  is  to  be  noted  that  the  {c^}  sequence  will  be  identically  zero  out¬ 
side  the  time  range  0  £  k  <_  s.  Upon  taking  the  discrete  Fourier 


transform  of  relationship  (4. 4. 1.4),  we  have  Sx(w)  in  the  form 


s>> 


co  +  « 


-jw 


.  +  c  e 
s 


-jsw 


1  +  a^  e-^U  + 


.  +  a  e 
P 


•jpm 


(4. 4. 1.5) 


If  this  expression  is  substituted  into  relationship  (4. 4. 1.3),  the 
required  formulation  of  the  spectral  density  estimate  is  completed. 


4.4.2  Smoothed  Periodgram  Method  (Cadzow,  1980  b) 

In  the  smoothed  periodgram  method,  one  first  generates  the 
auxiliary  "residual”  time  series  elements  according  to  the  relationship 

P 

e(k)  ■  x(k)  +  Z  a.  x(k  -i)  p  +  1  <_  k  £  n  (4. 4. 2.1) 

i-1  1 


in  which  the  ARMA  model’s  a^  coefficients  as  generated  by  relationship 


(4.3.7)  is  utilized.  Upon  examination  of  relationship  (4.2.1)  and 
under  the  condition  that  the  time  series  being  characterized  is  an 


ARMA  model  of  order  p  with  the  calculated  a^  coefficients,  it  follows  j 

that  the  residual  time  series  will  have  a  moving  average  spectral 
density  as  given  by 

S  (to)  =  J  l  b  e“jk“|2  (4. 4. 2. 2) 

k=0  I 

l 

i 

This  observation  in  conjunction  with  the  ARMA  model  representation 
then  provides  the  vehicle  for  estimating  the  underlying  time  series 
spectral  density,  that  is 

Sx(oj)  =  Se(cu)  /  |  l  ak  e"jktu|2  ,  aQ  =  1  (4. 4. 2. 3) 

k=0 

With  this  in  mind,  the  final  step  of  the  spectral  estimation  procedure 
requires  fitting  a  q-th  order  moving  average  (MA)  model  to  the 
residual  time  series  segment  (4. 4. 2.1)  to  effect  an  estimate  of  Se(w). 

St 

The  approach  to  be  presented  for  obtaining  the  q+1  order  MA 
model  is  an  adaption  of  the  well-known  method  of  Welch  for  obtaining 
smoothed  periodgrams  (Welch,  1967).  In  essence,  one  first  segments 
the  calculated  residual  elements  (4. 4. 2.1)  into  L  segments  each  of 
length  q+1  according  to 

et(k)  *  w(k)  e(k  +  1  +  p  +  id)  (4. 4. 2. 4) 

0  <  k  <_  q 
0  <  i  <  L  -  1 


v 
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where  w(n)  is  a  data  window  and  "d"  is  a  positive  integer  which 
specifies  the  time  shift  between  adjacent  segments.  These  individual 
segments  are  seen  to  overlap  for  a  shift  selection  of  d  q. 

Furthermore,  in  order  to  include  only  the  observed  time  series 
elements,  the  relevant  parameter  must  be  selected  so  that  p  +  q  + 

(L  -  l)d  <  n.  Finally,  the  q  +  1  order  periodogram  of  each  of  the 
L  segments  (4. 4. 2. 4)  is  taken,  and,  these  periodograms  are  in  turn 
averaged  to  obtain  the  desired  smoothed  q  +  1  order  MA  estimate  given 
by 

S  (w)  -  f  E  {-^T  I  2  w(k)  e(k  +  1  +  p  +  id)e_jwk  |2  }  (4. 4. 2. 5) 

e  L  .  „  q+1  ,  „  1 

i*0  M  k=0 

2 

where  the  data  window  is  normalized  according  to  E  w  (k)  *  1. 

In  using  this  smoothing  procedure,  the  variance  of  the  estimate 
Se(io)  is  reduced.  The  price  paid  for  this  reduction,  however,  is  a 
loss  in  frequency  resolution  and  an  increased  bias  of  the  estimate. 
Fortunately,  the  basic  resolution  capability  of  this  and  other  ARMA 
model  procedures  is  primarily  influenced  by  the  autoregressive  co¬ 
efficient  selection.  If  one  is  mainly  interested  in  resolution 
performance,  an  examination  of  the  ARMA  models'  pole  locations  then 
need  be  investigated. 

4.5  Numerical  Examples 

In  this  section,  the  classical  problem  of  detecting  the  presence 
of  sinusoids  in  additive  noise  is  considered.  In  particular,  we 
will  investigate  the  specific  case  in  which  the  time  series  observations 


are  generated  according  to 


x(n)  =  cos  (irf^n)  +  A2  cos  (irf2n)  +  w(n)  (4.5.1) 

1  _<  n  £  N 

where  w(n)  is  a  white  Gaussian  time  series  with  variance  one.  This 
particular  problem  serves  as  an  excellent  vehicle  for  measuring  a 
spectral  estimator's  performance  relative  to:  (i)  detecting  the 
presence  of  sinusoids  in  a  strong  noisy  background,  and  (ii)  resolving 
two  sinusoids  whose  frequencies  f^  and  are  nearly  equal.  The 
individual  sinusoidal  signal-to-noise  ratios  (SNR)  for  the  above  signal 
are  given  by  20  log  (A^//2)  for  k  *  1,2.  In  order  to  consider  the 
effectiveness  of  the  high  performance  ARMA  spectral  estimator  in 
different  noise  environments,  we  shall  consider  two  cases.  These 
cases  have  been  examined  in  reference  (Sullivan,  etc.,  1978)  where  the 
performance  of  many  modem  spectral  estimators  are  empirically  compared. 

CASE  I:  Ax  =  J20,  f  =  0.4 

A2  =■  /2 ,  f2  =  0.426 

In  this  example,  we  have  two  closely  spaced  (in  frequency)  sinu¬ 
soids  for  which  the  stronger  sinusoid  has  a  SNR  of  10  dB  while  the 
weaker  sinusoid  has  a  SNR  of  0  dB.  For  this  relatively  low  SNR  case, 
the  ability  of  a  spectral  estimator  to  resolve  closely  spaced  sinusoids 
and  identify  their  frequencies  will  be  tested.  Upon  generating 
sequence  (4.5.1)  with  the  postulated  parameters  for  a  data  length  of 
N  *  1024,  spectral  estimates  were  obtained  using  a  12-th  order  model 


with  the  high  performance  ARMA  method  (diagonal  element  of  the 
2 

matrix  A  is  (N-m)  ),  maximum  entropy  method,  and  the  Box-Jenkins 
method  incorporating  biased  autocorrelation  estimates.  In  addition, 
a  standard  periodgram  spectral  estimate  was  obtained  using  the  same 
data.  The  resultant  spectral  estimates  are  displayed  in  Fig.  4.5.1 
where  a  number  of  observations  can  be  made 

(i)  The  indirect  ARMA  spectral  estimate  provides  excellent 
results  with  two  sharp  peaks  at  f^  =  0.400  and 
f2  =  0.427,  and  with  the  spectrum  near  0  dB  (the  noise 
level)  for  most  other  frequencies. 

(ii)  The  maximum  entropy  and  Box-Jenkins  methods  were  unable 
to  resolve  the  two  sinusoids  in  the  prevailing  low 
SNR  environment. 

(iii)  Although  the  periodgram  is  able  to  resolve  the  two 

sinusoids,  the  well-known  random  fluctuation  behavior 
which  characterizes  the  periodgram  method  is  in 
evidence. 

This  example  nicely  demonstrates  the  potential  capability  of  the  high 
performance  ARMA  spectral  estimation  method  relative  to  existing 
procedures . 

In  many  practical  problems,  one  does  not  have  available  exceedingly 
long  data  lengths  upon  which  to  make  a  spectral  estimate.  To  demon¬ 
strate  the  ability  of  the  high  performance  ARMA  spectral  estimator  to 
perform  in  such  situations,  the  first  64  samples  of  the  data  sequence 
in  the  above  example  were  used  to  generate  a  spectral  estimate.  The 


Fig.  4.5.1(b)  Maximum  Entropy  Method  (P  “  15,  N  ■  1024) 


Modified  Box-Jenkins  Method  (P  =  15,  N  «*  1024) 
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resultant  15-th  order  high  performance  ARMA  spectral  estimate  obtained 
is  shown  in  Fig.  4.5.2  where  the  ability  to  resolve  the  two  closely 
spaced  sinusoids  is  again  evident.  The  sinusoid's  frequency  estimate 
f^  =  0.399  and  f9  =  0.423  are  also  of  good  quality  in  this  low  SNR 
environment . 


CASE  II:  Ax  =  /!,  f  =  0.32812 


A2  =  SZ,  f2  -  0.5 


We  are  now  examining  the  ability  of  the  ARMA  spectral  estimator 

to  detect  sinusoids  in  a  low  SNR  environment.  For  a  selection  of 

2 

N  =  64,  w(n)  =*  (N  -  n)  and  p  =  5,  the  resultant  ARMA  spectral 
estimation  is  displayed  in  Fig.  4.5.3(a).  Clearly,  one  is  able  to 
detect  the  presence  of  the  two  sinusoids,  and,  the  frequency  estimate 
f  =  0.3202  and  f2  =  0.5012  are  of  good  quality  considering  the 
prevailing  SNR  environment.  A  15-th  order  maximum  entropy  spectral 
estimator  was  then  found  to  generate  the  spectral  estimate  displayed 
in  Fig.  4.5.3(b).  Although  the  two  sinusoids  were  properly  detected, 
a  number  of  false  peaks  are  in  evidence. 

Next,  we  treat  the  time  series  recently  considered  by  Bruzzone 
and  Kaveh  (1980).  Specifically,  their  ARMA  time  series  is  characterized 
by 


where 

by 


Xk  =  \  +  °-5  £k 


1  2 

the  time  series  and 


(4.5.2a) 


are  autoregressive  process  generated 


Normalized  Frequency 


51 


x^  =  0.4  x^_x  -  0.93  x^_,  +  £*  (4.5.2b) 

xj^  »  -0.5  x^  ^  -  0.93  +  £k  (4.5.2c) 

1  2 

in  which  the  £,  ,  e,  and  e,  are  uncorrelated  Gaussian  random  variables 
k  k  k 

with  zero  mean  and  variance  1.  The  spectral  density  of  the  above 
time  series  (4.5.2a)  is  given  by 

Sx(oj)  =  1 1  -  0.4  e"ju)  +  0.93  e~j2olp2 

+  (l  +  0.5  e"j“  +  0.93  e_j2T2  +  0.25  (4. 5. 2d) 

Using  this  time  series  (4.5.2a),  twenty  different  independent  sampled 
sequences  each  of  length  64  were  generated.  These  twenty  observation 
sets  were  used  to  test  various  spectral  estimation  methods.  In 
Fig.  4.5.4,  twenty  superimposed  plots  of  the  ARMA  model  spectral 
estimates  of  order  (4.4)  as  obtained  by  using  the  Box-Jenkins  method,  the 
high  performance  method  with  t  5  4,  8  and  20  are  shown.  For  comparison 
purposes,  the  ideal  spectrum  is  also  plotted.  Comparing  the  two  cop 
most  plots,  the  high  performance  method  with  the  minimal  value  of  t = 4 
was  found  to  yield  a  marginally  better  spectral  estimate  than  the 
Box-Jenkins  method.  In  the  lower  two  plots,  one  can  observe  that  the 
high  performance  spectral  estimates  improve  significantly  as  t  is 
increased.  Next,  twenty  different  samples  sequence  of  length  200  were 
generated  according  to  time  series  (4.5.2a).  With  this  longer  data 
length,  it  was  anticipated  that  an  Improvement  in  spectral  estimation 
performance  would  result.  As  shown  in  Fig.  4.5.5,  a  marked  Improvement 


54 


is  evident,  where  the  ARMA  model  spectral  estimates  of  order  (4.4) 
are  shown  for  the  Box-Jenkins  method  and  the  high  performance  method 
for  selections  of  t  =  4,  8,  and  20. 

It  is  also  possible  to  use  the  high  performance  ARMA  method  for 
synthesizing  digital  filters.  To  illustrate  the  approach  that  is 
taken,  let  us  consider  the  specific  case  of  designing  a  low-pass 
filter  of  normalized  cutoff  frequency  f^.  One  may  readily  show  that 
the  impulse  response  of  an  idealized  version  of  this  low  pass  filter 
is  given  by  sin  (irfcn)/irn.  With  this  in  mind,  one  then  applies  the 
herein  developed  ARMA  procedure  to  the  specific  sequence 

x(n)  =  sin  [iTfc(n  -  0.5  N)]/ir(n  -  0.5  N)  1  <.  n  <_  N  (4.5.3) 

The  resultant  ARMA  model  obtained  in  this  manner  will  have  attenuation 
characteristics  of  the  desired  low-pass  filter.  To  illustrate  this, 
a  15-th  order  ARMA  spectral  estimate  of  this  sequence  was  made  for 
fc  =  0.2,  N  =»  128  and  w(n)  =  (N-n)  .  The  resultant  filter's  magnitude 
characteristics  are  displayed  in  Fig.  4.5.6  where  the  low-pass 
characteristics  are  In  evidence. 

4.6  Summary 

The  "high  performance"  ARMA  model  spectral  estimation  has  been 
described.  This  estimation  approach  provided  an  excellent  spectral 
estimation  performance  when  compared  with  such  contemporary  procedures 
as  the  maximum  entropy  and  Box-Jenkins  Methods.  The  above  mentioned 
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Performance  ARMA  Method  (P  =  15) 


Chapter  5 


COMPUTATIONALLY  EFFICIENT  ARMA  SPECTRAL  ESTIMATION 


5.1  Introduction 


Recently,  much  attention  has  been  focused  on  developing  spectral 

estimation  algorithms.  Unfortunately,  direct  application  of  the  linear 

prediction  method  as  described  in  Section  3.3.1  results  in  an 

excessive  computational  requirement,  since  it  is  necessary  to  solve 

3 

a  pxp  matrix  system  of  equations  which  generally  requires  0(p  ) 
computations.  For  this  reason,  a  number  of  computationally  fast 
algorithms  have  been  developed  to  overcome  this  difficulty.  These 
include  the  Levinson  algorithm  (Levinson,  1947).  The  Levinson 
algorithm  is  found  to  be  dependent  on  the  Toeplitz  structure  of 
the  matrix  characterizing  the  system  of  equations.  With  this  very 
restrictive  constraint  in  mind,  Kailath,  etc.  developed  the  concept 
of  the  displacement  rank  so  as  to  yield  efficient  solutions  for  non 
Toeplitz  system  of  equations.  The  displacement  rank  measures  how 
"close"  to  Toeplitz  a  given  square  matrix  is  (Kailath,  etc.,  1979). 

If  a  given  matrix  T  is  Toeplitz,  then  its  structure  is  characterized 
by  the  following  property 


T 


]  =  [t 


i+m, 


(5.1.1) 


where  t  denotes  the  (i,j)-th  element  of  the  pxp  Toeplitz  matrix 
1  J 

T  and  m  is  a  scalar  integer  (1  <_  i+m,  j+m  <_p).  That  is,  the  elements 
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of  the  matrix  T  are  identical  along  the  diagonal  and  subdiagonal 
directions.  In  recognition  of  this  key  property  of  a  Toeplitz 
structure,  the  displacement  rank  of  the  pxp  matrix  A  is  defined  by 


a (A)  =  min{a+  (A),  a  (A)} 

where 

a+  (A)  =  rank  {A  -  S  A  ST> 
a  (A)  =  rank  {A  -  ST  A  S} 


(5.1.2a) 


(5.1.2b) 


in  which  ct+  (A)  and  a_  (A)  are  called  the  positive  and  negative 
displacement  ranks  of  matrix  A,  respectively,  and  S  denotes  the  pxp 
down  shift  matrix  defined  by 


(5.1.3) 


It  can  be  straightforwardly  shown  that  the  displacement  rank  of  a 
Toeplitz  matrix  T  is  2,  that  is 


a(T)  =  a+  (T)  =  (T)  =  2  (5.1.4) 

If  a  given  matrix  A  has  a  displacement  rank  a,  then  it  has  been 

shown  that  the  Inversion  of  A  may  be  accomplished  with  the  number  of 

2 

required  computations  being  0(ap  )  (Friedlander ,  etc.,  1979). 

Based  on  these  concepts,  a  number  of  computationally  efficient 
algorithms  for  AR  spectral  models  have  been  developed  (Friedlander, 
etc.,  1978,  1979;  Morf,  etc.,  1977;  Morf  and  Lee,  1978;  Lee  and  Morf, 
1980;  Morf  and  Kailath,  1975;  Mullis  and  Roberts,  1976;  Morf,  1980; 
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Bimead  and  Anderson,  1979)  .  Some  of  these  methods  are  classified 
by  Morf,  etc.  (Morf,  etc.,  1977). 

In  this  chapter,  fast  algorithms  which  are  applicable  to  the 
"high  performance”  ARMA  method  (see  Chapter  4)  are  developed.  To 
achieve  the  fast  algorithm  solution  capability,  it  will  be  necessary 
to  restrict  the  number  of  Yule-Walker  approximation  to  be  p  (i.e., 
t  =  p) .  Unfortunately,  the  restriction  t  =  p  will  generally  result 
in  an  associated  decrease  in  spectral  estimation  performance.  Thus, 
in  obtaining  a  computationally  fast  algorithmic  solution  procedure 
for  the  a^  coefficients,  an  accompanying  sacrifice  in  spectral 
estimation  performance  is  the  price  being  paid.  One  must  therefore 
carefully  consider  the  tradeoff  for  any  given  application.  Fortunately, 
the  degradation  in  performance  is  not  great  for  many  relevant 
applications  in  which  the  data  length  n  adequately  exceeds  the  ARMA 
model  order  parameters  p  and  q. 

The  achievement  of  fast  algorithms  requires  data  modifications 
which  will  be  discussed  in  Section  5.2.  In  Sections  5.3  and  5.4, 

2 

algorithms  which  requires  0(p  )  and  0(p  log  p)  multiplications, 
respectively  are  discussed.  An  algorithm  which  requires  0(p) 
computations  is  developed  in  Chapter  6. 

5.2  Data  Modification 

In  this  section,  we  will  discuss  three  types  of  data  modifications 
referred  to  as  the  pre-modification,  post-modification  and  pre-  and 


6' 


W 


post-modification  methods  (Cadzow  and  Ogino,  1981) .  These  are  modi¬ 
fications  of  the  'high  performance'  ARMA  spectral  estimation  methods 
as  discussed  in  Section  4.3  in  which  t  is  restricted  to  be  p.  It  will 
be  recalled  that  in  this  unmodified  case  one  must  solve  the  matrix 
system  of  equations  (4.3.4).  Without  loss  of  generality,  this 
matrix  system  of  equations  may  be  represented  as 

Y+  X  a  =  Y*  x  (5.2.1) 


where  Y  and  X  are  (n  -  p)  x  p  Toeplitz  matrices,  while  x  and  a_  are 
(n  -  p)  x  1  and  p  x  1  column  vectors,  respectively  defined  by 


t  >  X  .  -  , 

p-q  p-q+i 


l-q  2-q’ 


t 


)  A  1 

n-q-1 


n-q-p 


(5.2.1a) 


x  x  ,, 

P  p+1 


X1  x2 


.  .  x 


h-1 


n-p 


£  "  [xp+l’  xp+2’  ’  •  *  ’  xn] 


(5.2.1b) 

(5.2.1c) 


=  Ca. ,  .  •  .  ,  a  1 
1  P 


(5. 2. Id) 


where  the  entries  of  the  matrices  X,  Y  and  column  vector  x  can  be 
determined  from  expression  (4.3.4).  The  entries  of  the  column  vector 
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a  in  expression  (5.2. Id)  denote  the  denominator  AR  coefficients  to  be 

found.  It  can  be  shown  that  the  displacement  rank  of  the  matrix 

YTX  is  4.  As  suggested  in  Section  5.1,  it  is  possible  to  find  an 

2 

algorithmic  solution  procedure  which  requires  0(4p  )  computations. 

In  fact,  in  Section  5.3,  a  generalized  Levinson's  algorithm  will  be 
developed. 

It  is  possible  to  realize  significant  computational  savings  in 
the  'high  performance'  ARMA  spectral  estimation  procedure.  This 
improvement  will  entail  a  slight  modification  in  the  vector  x  and 
matrices  X  and  Y.  Although  the  suggested  modifications  will  typically 
result  in  biased  estimates  of  the  Yule-Walker  equations,  it  is  shown 
that  when  the  data  length  n  adequately  exceeds  the  order  parameter 
p  and  q  then  these  estimates  are  virtually  unbiased  (Cadzow  and  Ogino, 
1981) . 

With  the  above  high  performance  spectral  estimation  method 
representation  serving  as  a  basis,  we  shall  now  consider  the  afore¬ 
mentioned  modifications  required  to  achieve  computationally  efficient 
algorithmic  solution  procedures. 

5.2.1  Pre-modification  Method 

In  expressions  (5.2.1a)  and  (5.2.1b),  the  addition  of  lower 
triangular  matrices  to  the  top  of  matrices  X  and  Y  yields  the  Toeplitz 


matrices 
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n-q-1 


n-q-p 


(5.2.1.1a) 


T 


n-1 


n-p 


(5.2.1.1b) 


with  Y^  and  X^  each  being  (n  -  1)  x  p  matrices.  While  maintaining 
the  structure  of  expression  (4.3.1),  the  vector  x  will  be  modified 
to 


C^c-j  *  •  •  (5.2.1.1c) 

Substitution  of  expressions  Y^,  and  x-^  in  place  of  Y,  X  and  x, 
respectively,  yields 

Y^X1a=Y^x1  (5.2.2) 

It  can  be  shown  that  the  displacement  rank  of  the  matrix  Y^  X^  is  3 

(Cadzow  and  Ogino,  1981).  It  is  possible  to  find  a  generalized 

2  t 

Levinson  algorithm  which  requires  0(3p  )  computations  to  invert  Y^X^. 
More  importantly,  because  of  this  specific  structure,  an  algorithm 
which  requires  0(p)  computations  has  been  developed  and  will  be 
discussed  in  Chapter  6. 


5.2.2  Post-modification  Method 


Following  a  similar  procedure  as  employed  in  Section  5.2.1, 
the  addition  of  an  upper  triangular  matrix  to  the  main  body  of  the 
matrices  specified  by  (5.2.1. a)  and  (5.2.1b)  yields  the  Toeplitz 
matrices 


*  « 
n-q-1 


•  o 


X  ....  X 

n-q-p  n-q-1 


T 


(5.2.1.3a) 


where  and  are  each  (n  -  1)  x  p  Toeplitz  matrices.  In  a  similar 


manner,  the  column  vector  x^  is  defined  by 


(5.2.1.3c) 


The  displacement  rank  of  the  matrix  Y^  ^  is  f°und  to  be  (Cadzow 
and  Ogino,  1981) 


a+  (Y2  V  *  a-  (Y2  X2)  *  3 


(5. 2. 1.4) 


r 
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thereby  offering  a  generalized  Levinson  solution  procedure  requiring 

2 

a  computational  complexity  of  0(3p  )  for  solving  the  system  of 
equations 


Y2  X2  a 


Y2  x 


(5. 2. 1.5) 


A  more  computationally  efficient  algorithm  associated  with  the  post¬ 
modification  will  be  developed  herein.  It  is  shown  that  the  number 
of  computations  is  reduced  to  0(p  log  p)  if  p  =  q  where  p  and  q  are 
the  order  of  denominator  and  numerator  coefficients  of  the  ARMA  model, 
respectively. 

5.2.3  Pre-  and  Post-Modification  Method 

The  combination  of  the  previously  discussed  two  modification 
methods  yields  the  pre-  and  post-modification  method.  The  matrices 
and  vector  are  modified  in  the  following  manner 


x-  •  •  •  x  •  • 
1-q  p-q 


I  T 


*  » 
n-l-q 


o 


C\  •• 

W  xl-q  *  *  *  xn-p-q 


(V  - 

n-l-q 


(5.2.3.1a) 


•  *  x  *  • 

P 


."-l  .  o 

•  »  •  • 

o 

X.  ...  X  ...  .  x  , 

l  n-p  n-1 


(5.2.3.1b) 


65 


x-  =  [x«  .  .  .  x  ,  ...  x  0 . 0]  (5.2.3.1c) 

— J  ^  p+i  n  s _ _ _ 

p 

rzeros 

where  and  denote  (n  +  p  -  1)  x  p  matrices  and  is  a 
(n  +  p  -  1)  x  1  column  vector  respectively. 

It  can  be  shown  that  the  matrix  X^  is  a  Toeplitz  matrix.  A 
conventional  approach  for  solving  the  Toeplitz  system  of  equations 

Y3  X3  a  =  Y^  x3  (5.2.4) 

2 

was  developed  by  Levinson  (Levinson,  1947),  which  requires  0(p  ) 
computations.  More  recently  much  effort  has  been  conducted  in 
developing  more  efficient  AR  algorithms  whose  computational  require¬ 
ment  is  0 (p  log  p) .  Gustavson,  etc.,  presented  their  algorithms 
which  were  based  on  the  use  of  Pade  approximates  and  the  rational 
Hermite  approximation  (Gustavson  and  Yun,  1979).  Morf 
developed  the  so-called  doubling  algorithm  which  requires 
0(p  log  p)  (Morf,  1980).  Bitmead  and  Anderson  also  independently 
found  a  doubling  algorithm  (Bitmead  and  Anderson,  1979).  In  Section 
5.4,  an  application  of  the  doubling  algorithm  to  the  ARMA  model  is 
developed. 

5.3  Generalized  Levinson's  Approach  for  the  ARMA  Model: 

The  Unmodified  Method 
» 

In  this  section,  an  algorithm  which  can  be  applied  to  the  direct 

approach  (i.e.,  no  modification)  will  be  developed.  Without  loss  of 

generality,  the  m  x  m  matrix  R,  will  be  defined  by 

1 
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r“  -  [y"  ]t  X? 

1,01  l,m  l,m 


where  X,  denotes  (n  -  m  +  1)  x  m  matrix  defined  by 
l,m 


(5.3.1a) 


(5.3.1b) 


with  the  subscript  m  designating  the  number  of  columns  of  matrix 

X?  ,  1  is  the  smallest  and  n  the  largest  index  of  the  observation 
l,m 

data  to  form  the  matrix  X?  .  In  a  similar  manner,  matrix  y”  is 

l,m  l,m 


obtained  by 


where  the  entries  of  the  matrix  y”  are  given  by 

l,m 

it 

yi  m  for  i  *  If  •••  ,  n 


(5.3.1c) 


(5. 3. Id) 


This  particular  representation  has  been  chosen  so  that  in  the  develop¬ 
ment  of  the  generalized  Levinson's  algorithm  for  an  ARMA  model, 
notational  complexity  can  be  eased.  It  then  follows  that  the  matrix 
expressed  by  (5.3.1a)  has  the  following  shift  invariance  structures 

which  characterizes  the  near  Toeplitz  structure  of  the  matrix  r!?  , 

l,m 


that  is 
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_n  „n-l  ,  n  .  n.T 

R,  =  R,  +  y_  (x  ) 
l,m  l,m  -th  — m 


(5.3.2a) 


_n  ,  m  ,  m.T 

R,  +  v  (x  ) 

2,111  — m 


(5.3.2b) 


l.mfl 


,  n  .T 

(W 


<il)T 


(5.3.2c) 


(5. 3. 2d) 


where  x^»  ^  and  x^  are  m  x  1  column  vectors  defined  by 


i  =  [V  •  *  *  ’  yn-nH-l]1 


(5 . 3. 2e) 


nr  -i" 

x  =Lx,...,x  ,,J 

-ra  n  n-m+1 


(5.3.2f) 


ra  r  -il 

2m  *  CV  *  •  *  ’  yl] 


(5.3.2g) 


m  r  il 

5m  m  LX  t  .  .  .  >  X-  J 
-ra  m  1 


(5.3.2h) 


while  (z^^)T  and  (w^^)^  denote  the  first  and  last  rows  of  the 

matrix  R^  -  respectively,  and  the  m  x  m  matrices  R?  is  defined  by 
X  f  ot  l  z 


_  -  [y”  m]T  x" 

z  ,m  z  ,m  z  ,m 


(5.3.3) 
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From  a  structural  viewpoint,  relationship  (5.3.2a)  is  called  time 


In  expression  (5.3.2c),  R^  ^  is  seen  to  be  a  submatrix  of  r” 

l,m  l,m+l 


update,  since  the  matrix  R,  is  explicitly  defined  as  a  sum  of  two 

l,m 

matrices: R^  ^  which  includes  all  of  the  past  observations  up  to 
previous  time  index  n-1  and  y^(x^)^  which  includes  the  most  recent  data. 

It 

then  follows  that  relationship  (5.3.2c)  is  called  order  and  time 
update. 

A  computationally  efficient  algorithm  will  be  obtained  by  using 
various  combinations  of  the  above  shift  invariance  structures.  This 
fast  algorithmic  procedure  for  finding  the  solution  is  similar  to 
Levinson's  algorithm  (Levinson,  1947).  The  overall  solution  is 
updated  from  the  solution  of  a  lower  order  to  that  of  higher  order 
system  of  equations  (order  update)  and  from  the  solution  of  previous 
time  instance  to  that  of  present  time  (time  update).  To  develop  this 
algorithm,  we  apply  an  induction  hypothesis.  Suppose  at  order  m 
and  time  n,  we  have  the  relationship 


_n  r  n  .  n  ,n  n 

R.  La-,  b,  d,  I  = 

l,m  l,m  — l,m  — l,nr 


e,n 

n,m 


m 


.r,n 


(5.3.4a) 


J  l»m  '1| 

where  a“  ,  b“  ,  and  d“  are  m  x  1  column  vectors  defined  by 

1  9  ID  X  9  ttl  1  9  01 

a"  *  [1,  a?  (1),  .  .  .  ,  a"  (m-l)]T  (5.3.4b) 

X9II1  X9Q1 


b"  -  [b?  (n-1) ,  b"  (n-2) , 

X  9  IQ  X  9  01  X  9  01 


n 

>l,m' 


bV  (1),  l]1  (5.3.4c) 
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and 


d“  -  [d“  (1),  d“  (2) . d"  (m)]1 

— l,m  l,m  1,111  l,m 


(5.3.4d) 


Specifically,  a  and  b.  are  called  the  forward  and  the  backward 
— 1 ,  m  — 1 ,  m 

AR  coefficient  vectors,  respectively.  In  the  development  of  the 

computationally  efficient  algorithm,  the  auxiliary  vectors  d”  are 

— I,® 

needed  to  cancel  the  end  effects  due  to  the  non-Toeplitz  nature  of 


matrix  R  .  At  the  previous  time  index  n-1,  we  have  the  relationship 
x  f  m 


nn-l  r  n-1  .n-1  .n-li 

Ri  Lai  b.  d  J 

l,m  l,m  — l,m  —  l,m 


.e,n-l 


’l,m 


m 


„r,n-l 


"l,m 


(5.3.4e) 


Based  on  the  relationship  (5.3.4a)  and  (5.3.4e),  we  will  develop  a  recur- 

n  ,  n  ,  .n 


sive  solution  procedure  for  the  vectors  a,  . , ,  b,  . .  and  d.,  . , 

J. )  IQ *  X  X )  m  *  1  X  j  in  I  1 

function  of  n.  Applying  the  shift  invariance  structure  (5.3.2b)  to 
(5.3.4a)  yields  the  following  expressions 


as  a 


R?  a? 
2,m  — l,m 


e,n  m 

5.  e  -  e  y_ 
l,m  — 1  m  “Tn 


(5.3.5a) 


i,  d“  -  (1  -  f  )  y 
2,m  — l,m  m  -Tn 


(5.3.5b) 


where  e  and  f  are  scalars  defined  by 
mm  J 


,  m.  T  n 
e  -  (x  )  a, 
in  “in  x  f  in 


(5.3.6a) 


/  m\T  .n 

(x )  d 

—m  — l,m 


(5.3.6b) 
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and  &  is  the  m  x  1  unit  basis  vector  expressed  by 

e^  “  [l,  0  .  .  .  o]T 


Expressions  (5.3.5a)  and  (5.3.5b)  lead  to  the  following 
relationship 


n  ,.n 

R2,m  — 2,m  "  C2,m  -1 


(5.3.6c) 


(5.3.7a) 


where  a°  and  are  a  m  x  1  column  vector  and  a  scalar  respectively, 

6  j  m  4  j  m 


defined  by 


4  =  C4  + ,  4  m  i/{1  +  -f-'V  4  J1)}  <5-3-7b> 

^  j  in  l,m  X  “  t  1  j  in  l""*!  l,m 


5, ’  *  c,*  /  (l  +  i - r-  d.  m(i)} 

l,m  1  -  t  l.m 

m  * 


(5.3.7c) 


where  d^  (1)  denotes  the  first  entry  of  the  vector  d”  (see 

1  f  ®  l,m 

eq.  (5.3. 4d)).  Expressions  (5.3.2c)  and  (5. 3. 2d)  lead  to  the  relation¬ 
ships 


_n  n  e,n 

l.nri-l  — l,nrfl  “  5l,nri-l  -1 


(5.3.8a) 


„n  ,n  =  r,n 
l.nrfl  l.mfl  =  ?l,nrt-l  ^nri-l 


(5.3.8b) 


where  e  . ,  is  the  (nrfl)xl  unit  bases  vector  defined  by 
— nrrl 


Co,  ....  i]1 


(5.3.8c) 


In  expressions  (5.3.8a)  and  (5.3.8b),  ,n. ,  and  c5,n. ,  are  scalars 

i,nn-i  i,nrri 

defined  by 
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e,n  _ 
^l,mfl  z,m 


a  8 
m  m 

r,n-l 

^l,m 


(5.3.9a) 


i  a  B 

r,n  m  r,n-l  _  m  m 
5l,nrt-l  5l,m  c,n 

52,m 


in  which  a  and  8  are  scalars  specified  by 
m  m 


,  n  ,T 

“m  *  <W 


n 

— 2,m 


(5.3.9b) 


(5.3.10a) 


m 


/  n  nT 


b“_1 

— l,m 


(5.3.10b) 


The  m  x  1  column  vectors  a?  and  b"  , ,  in  expressions  (5.3.8a) 

—l,mrl  — l,nrri 

and  (5.3.8b),  respectively,  are  defined  by 


■l.nrt-l 


- 

a 

-  - 

m 

m 

0 

— 2,m 

o 

r,n-l 

n,m 

»r1 

— l,m 

m 

(5.3.11a) 


hn 

— l,m+-l 


0 

8 

•  ■ 

m 

n 

b?_1 
— l,m 

e,n 

^2,m 

— 2,m 

0 

ft 

■  ■ 

(5.3.11b) 


Expressions  (5.3.11a)  and  (5.3.11b)  are  seen  to  be  very  similar  to 

Levinson's  algorithm  (Levinson,  1947).  In  fact,  one  can  show  that 

these  two  expressions  can  be  converted  to  Levinson's  algorithm,  if  the 

pre-  and  post-modification  method  is  applied  on  the  matrix  r!?  . 

Next,  we  will  verify  the  relationship  which  updates  the  vector 

d?  .  The  m  x  1  column  vector  d^  . ,  is  found  to  be 
— l,m  — l.ntfl 
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jn 

-l,m+l 


0 

d“ 

— l,m 


ym+l  '  Ym  n 


,e,n 

’l,m+l 


— l.nrfl 


(5.3.12a) 


where  y  is  a  scalar  given  by 
m 


,  n  .T 

Ym  *  <W 


d“ 

— l,m 


It  can  be  straightforwardly  shown  that 


(5.3.12b) 


_n  ,n  _  m+1 

±,mfl  -l,m+l  =  ^-m+1 


(5.3.13) 


Finally,  combining  expressions  (5.3.11a),  (5.3.11b),  and  (5.3.13), 
the  following  relationship  is  obtained. 


Ri,m+1  t-l.nH-1  -l.nH-1  -l.nrfl-* 


rE,n  0 

4l,nH-l 


0 

r 


Vfl 


0  Cl,m+1  yl 


(5.3.14) 


In  the  above  development,  the  generalized  Levinson's  algorithm  for 

AKMA  model  is  verified  based  on  the  induction  hypothesis.  The  number 

2 

of  computation  of  the  algorithm  is  readily  found  to  be  0(3p  )  where 

P  designates  the  number  of  denominator  coefficients  of  the  ARMA  model. 

We  will  now  detail  steps  of  the  computations  required  in 

this  recursive  algorithm.  The  algorithm  starts  with  the  initialization 

N 

procedure  at  n  *  q+1  and  order  m  ■  1.  The  solution  a 

— p+1 


\ 


of  the  matrix 
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equation  (5.3.4a)  with  m  =  p+1  is  obtained  by  recursively  updating 

a11  from  m  =  1  to  p  +  1  (order  update)  and  from  n  =  q+1  to  N  (time 
— m 

update).  Meanwhile,  auxiliary  vector  dn  is  also  recursively  updated. 

— m 

The  above  algorithm  can  be  presented  as  follows 


Step  1:  Initialization  for  time  update  (n  =  q+1) 


r  1  -i  r  1  -iT  r,l  _E»1 

51  ki.iJ  =  =  h,i  =  y 


q+l  Xq+1 


1  .  1  -  ,1  ;„£,1 

—1,1  —1,1  *  1  *  dl,l  ~  ^1^1,1 


Step  2:  n  M  n+1 


Step  3:  Initialization  for  order  update 


wn  =  wn^'  +  y  x 
— m  — m  n  n+l-m 


n  n-1  . 

z  =  z  +y^,wX., 

— m  — m  n+I-M  n+l-m 


where  M  =  min  (p+1,  n-q) 


m  =  1, 

m  =  1, 

m  *  1, 


,  M 

.  ,  M 

.  ,  M 


c *n  „r n  / .  v 
C, ’=?,’=  z  (1) 
i,m  i,m  — m 


where  zn(l)  denotes  the  first  element  of  column  vector  z11 
— m  — m 


n  ,  n  _  ,  .it  ,_e.,u 

1,1  ‘  bl,l  ■  1  di,i  ■  H'h.m 


e  ,n 


St0?  4:  Compute  recursively  from  m  =  1  to  M  where  M  -  min 
(p,  n-q-1) 


,  m.T  n  ,  ,  m.T  ,n 

£  m  (x  )  a,  ,  f  *  (x  )  d, 
m  — m  — l,m  m  -m  —  l,m 
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£  -  =  U?  _  +  d?  _]/{!  +  d?  _(!)} 


— 2,m  — l,m  1  -  f  — 1  m 

m 


1  -  f  -l,m' 
m 


/U  +  — T'  d?  (1)} 


’2,m  l,m 


1  -  f  — l,m 
m 


•  — 

* 

,  n  s  T 
a  =  (i..) 
m  — nrt-1 

-2,m 

.  o  _ 

*m  " 

0 

.r1 

— l,m 

■  m 

Update  forward  and  backward  solutions 


n 

-l.nri-l 


— l,nH-l 


n 

— 2,m 

0 


0 

n-1 


a 


.r,n-l 


l,m 

6 

m 

£.n 

Q2,m 


0 

b”-1 
— l,m 


n 

a_ 

-2,n 


a  6 

,e ,n  m  m 


e.n  _ _ 

l’2,m  r,n-l 

5l,m 

,  a  6 

r,n  _  r,n-l  _  m  m 
^l,m+l  ^l,ra  e,n 

?2,m 


Compute  auxiliary  vector  d™ 


m 


Ym  "  <^1>T 


d" 

— l,m 


Hn 

—1,13+1 


0 

d? 

— l,m 


ym+l  -  Ym 
e,n 

n,m+l 


n 

— l,nrfl 


75 


Step  5:  If  n  <  N  go  to  Step  2 
Step  6:  End  of  algorithm 

In  above,  N  is  taken  to  be  the  last  index  of  vector 


V 


5.4  ABMA  Doubling  Algorithm:  The  Pre-  and  Post-Method 

As  described  in  section  5.2,  one  of  the  data  modification 
methods  referred  to  as  the  pre-  and  post-modification  method  leads 
to  the  following  set  of  equations 


A  a  =  b  (5.4.1) 

p  -  ~p 


where  A  is  a  pxp  Toeplitz  matrix  and  b  is  a  pxl  column  vector 
P  ~P 

given  by 


A 

P 


T 

Y3  X3 


(5.4.2) 


b 

-P 


Yt  x 
3  -3 


where  matrices  Y^,  and  column  vector  x^  are  previously  defined  in 

expressions  (5.2,3.1a),  (5.2.3.1b)  and  (5.2.3.1c),  respectively.  The 

displacement  rank  of  the  matrix  A^  is  readily  shown  to  be  2.  This 

being  the  case,  it  is  possible  to  apply  the  doubling  algorithm  (Morf, 

1980;  Bitmead  and  Anderson,  1979). 

Without  loss  of  generality,  we  now  assume  that  p  *  2  for  some 

integer  k.  The  matrix  A  can  be  partitioned  into  4  matrices  whose 

P 

sizes  are  2^  ^  x  2^  ^ .  Each  2^  ^  x  2^  ^  matrix  is  then  also 
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(k-2)  (k-2) 

partitioned  into  2  x  2  matrices.  This  procedure  is  called 

the  doubling  or  halving  procedure.  In  this  procedure,  we  can  express 
a  21  x  21  matrix  A2^  in  terms  of  l  x  Jt  submatrices  B^,  C£,  and 

Ej  in  the  following  manner 

(5.4.4) 

and  its  inverse  is  found  to  be  the  form 


21 


(5.4.5) 


where  S^,  T^,  and  are  i  x  £  square  matrices  given  by 


s*  ’  h1  +  h1  c*  \  \  h1 


(5.4.6a) 


T.  =  -S„  C„  ET1 
K.  ill 


(5.4.6b) 


•  -E^  st 


v*  ■  EI1  +  'I1  D«  st  '«  'I1 


(5.4.6c) 

(5.4.6d) 


Relationships  (5.4.6a)  -  (5.4.6d)  are  straightforwardly  derived  from 
the  Schur  complements  theorem  (Aho,  etc.,  1974).  From  the  above 
relationships,  we  can  obtain  A^  from  and  E,1.  The  solution  of 
the  equation  (5.4.1)  requires  0(2  c(m))  computations  where  c(m)  is 
the  number  of  operations  required  to  multiply  a  vector  times  a 
triangular  Toeplitz  matrix. 
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The  number  of  computation  c(m)  is  obtained  in  a  following  manner. 
By  definition  (Kailath,  etc. ,  1979) ,  can  be  decomposed  in  the  form 


,£,  4  ”2* 

1=1 


(5.4.7) 


where  and  are  lower  and  upper  triangular  Toeplitz  matrices, 
respectively,  which  can  be  obtained  recursively  (Bitmead  and  B. 
Anderson,  1979) .  The  matrices  and  are  expressed  by 


i-ja.i) 

o 

Lj(2,l) 

Lj(2,2) 

uja.i) 

V\(1'2) 

o 

uj(2,2) 

(5.4.8a) 


(5.4.8b) 


where  L*(l,l)  and  L*(2,2)  are  l  x  l  lower  triangular  Toeplitz  matrices, 

U^(l,l)  and  U^(2,2)  are  l  x  l  upper  triangular  Toeplitz  matrices, 

and  L*(2,l)  and  U^(l,2)  are  l  x  i  full  Toeplitz  matrices.  Substitution 

of  expressions  (5.4.8a)  and  (5.4.8b)  into  (5.4.7)  yields  the 

partitions  of  the  matrix  A in  expressions  (5.4.4)  to  be 
2 

=»  ill  Lj(lfl)  uj(l,l)  (5.4.9a) 

2  2 

C£  =  lJ(1,1)  CuJ(1,2)]l  +  E^  lJ(1,1)  [^(l,2)]v 


(5.4.9b) 


2  .  2 
D*  -  2  ClJ(2,1)]  Uj(l,l)  +  E  [lJ(2,1)]  uj(l,l) 
i=l  i=l 


(5.4.9c) 
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2  2 

Et  =  2  l£(2,1)  uj(l,2)  +  2  lJ(2,2)  U*(2,2) 
i-1  i=l 


(5.4.9d) 


where  Che  following  relationships  are  implicitly  used 


UJ(1,2)  •  Oj«,2)]L  +  [0j(l,2)]0 


(5.4.10a) 


Lj(2.1)  =  [lJ(2,1)]l  +  ClJ(2,1)]u 


(5.4.10b) 


in  which  [u*(l,2)]  and  [L^(2,l)]^  denote  i  x  l  lower  triangular 
matrices,  and  [U*(l,2)]  and  [1^(2,!)]^  denote  i  x  l  upper  triangular 
matrices.  The  partitions  given  by  equations  (5.4.93)  -  (5.4.9d) 
are  expressed  in  terms  of  lower  triangular  and  upper  triangular 
matrices.  It  turns  out  that  the  use  of  above  relationships  reduces 
the  computational  complexity  c(m)  to  be  0(p  log  p)  (Morf,  1980). 

The  algorithm  which  makes  use  of  the  doubling  method  can  be  found  in 
(Morf,  1980).  Morf  described  the  algorithm  by  introducing  high 
computer  language  which  necessitates  frequent  subroutine  calls.  On 
the  other  hand,  the  step-wise  description  of  the  halving  method  is 
presented  in  (Bitmead  and  Anderson,  1979).  Implementation  of  the  halv¬ 
ing  algorithm  is  relatively  complex  and  a  rather  large  value  of  p  is 
required  before  the  computationally  complexity  0(p  log  p)  is 
approached. 


5.5  Numerical  Example 


In  this  section,  the  spectral  performance  of  the  pre-  and  the 
post-modified  methods  are  compared  with  the  unmodified  method.  As 
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a  test  example,  /e  treat  the  time  series  (4.5.2a).  Using  this  time 
series  (4.5.2a),  twenty  different  independent  sampled  sequences  each 
of  length  64  were  generated. 

The  modification  methods  were  applied  to  these  twenty  different 
sampled  sequences  of  length  64  to  obtain  ARMA  model  spectral  esti¬ 
mates  or  order  (4,4).  The  resultant  spectra  are  shown  in  Fig.  5.5.1. 
It  is  apparent  that  only  a  small  degradation  in  spectral  estimation 
performance  has  been  shown  by  the  modified  method.  It  might  be 
conjectured  that  the  implementation  of  the  fast  algorithms  will  not 
much  degrade  spectral  performance  in  many  practical  examples. 


Summary 


In  this  chapter,  computationally  efficient  ARMA  spectral 

estimation  algorithms  have  been  developed.  These  algorithms  are 

predicated  on  the  utilization  of  data  modification  methods. 

Specifically,  two  algorithms  referred  as  the  generalized  Levinson's 

algorithm  and  the  doubling  algorithm  were  developed  for  obtaining 

AR  coefficients  of  ARMA  model.  These  algorithms  have  a  computational 
2 

complexity  of  0(p  )  and  0(p  log  p),  respectively. 


Chapter  6 


A  RECURSIVE  ARMA  SPECTRAL  ESTIMATOR: 
THE  PREMODIFIED  METHOD 


6.1  Introduction 

A  recursive  ARMA  spectral  estimation  procedure  is  developed  in 
this  section.  It  is  recursive  in  the  sense  that  as  a  new  element  of 
the  time  series  is  observed,  the  parameters  of  a  spectral  estimation 
model  are  algorithmically  updated.  The  recursive  algorithm  requires 
0(p)  computations  to  update  the  model's  parameters  for  each  new  data 
point.  The  development  of  this  algorithm  is  predicated  on  utilization 
of  certain  projection  operators.  In  Section  (6.2),  a  vector  space  is 
formulated  by  making  use  of  the  given  observation  data.  The  method 
of  linear  predictions  will  give  rise  to  projection  operators  which 
decompose  relevant  vector  spaces  into  subspaces  spanned  by  the 
prediction  error  vector  and  the  observation  vectors.  Linear  prediction 
methods  used  in  this  chapter  include  forward  prediction,  backward 
prediction  and  delayed  backward  prediction.  Each  of  these  methods  is 
associated  with  its  own  projection  operator.  The  decomposition  of 
these  projection  operators  is  discussed  in  Section  (6.4).  The  order 
update  and  time  update  recursions,  as  described  in  Sections  (6.5)  and 
(6.6)  play  a  central  role  in  the  overall  recursive  algorithm.  Finally 
in  Section  (6.7),  a  recursive  algorithm  is  outlined. 


31 
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6.2  Vector  Space  Formulation 

In  this  section,  the  given  spectral  estimation  problem  will  be 
cast  into  a  convenient  vector  space  setting.  It  will  be  assumed  that 
the  following  observations  of  the  time  series  (x(n)} 


x!»  x2’ 

are  given.  This  in  turn  will  give  rise  to  the  associated  column 
data  vector 


(6.2.1) 


-  [x1,  x 2*  ...  ,  x^] 


(6.2.2) 


It  is  convenient  to  form  an  auxiliary  column  vector  specified  by 


%- s  % 


-  [0  ...  0  ... 


(6.2.3a) 

(6.2.3b) 


where  S  denotes  the  NxN  down  shift  matrix  (see  eq.  (5.1.3))  and q is  the 
numerator  order  of  the  ARMA  model.  The  vectors  x^  and  y^  lie  in  the 
product  space 


„N 

R  *  R  x  R  x  . . .  x  R 


(6.2.4) 


We  next  construct  the  subspace  which  is  spanned  by  the  set  of 

vectors  S*  x^»  •••  >  Sm  x^.  This  subspace  will  be  suggestively  denoted 

by 

M  Ml."]  ’  (S±  V  Sl+1  % . S“  V  <6'2-5) 
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where  the  first  integer  i  may  take  on  any  value  in  the  set  {0,  1,  ...  , 

m}.  As  will  be  described  in  Section  (6.3),  the  recursive  algorithm 

is  derived  for  particular  selections  of  indices  i  and  m.  Similarly, 

N 

for  the  vector  ^  contained  in  the  product  space  R  ,  the  associated 
subspace  M  ^  is  defined  by 

M  %[!,»]■ (s±  v  sl+1  % . S*V  <6-2-6> 

where  the  first  integer  i  may  take  on  any  value  in  the  set  {0,  1,  ...  , 
m}.  Next,  we  let  P  x^j-^  designate  the  projection  operator  on  the 
subspace  M  2^^  along  the  subspace  orthogonal  to  M  (this 

orthogonal  subspace  will  be  denoted  by  H  m]^'  ^is  projection 

operator  which  depends  on  x^  and  can  be  shown  to  have  the  form 

P  %[i,m]  =  A  %[i,m]  -%[i,m]  A  ^N[i,m]^  A  ^N[i,m] 

(6.2.7) 

where  A  x^j-^  m]  and  A  ^j[i  m]  are  ^  x(m-i+l)  matrices  composed 
of  the  following  ordered  set  of  column  vectors 

A  ’  Csl  %  sl+1%  •••  S”%]  <6-2'S) 

*%[!..]•  Csl  %  Sl+1%  •••  S"%]  (6-2'9) 

The  projection  characteristics  of  operator  (6.2.7)  are  depicted 
in  Fig.  6.2.1.  It  will  be  convenient  to  introduce  a  projection  opera¬ 
tor  on  the  complement  of  subspace  M  25^[^  m]*  This  operator  is  defined 
by 
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P  "  1  "  P 


(6.2.10) 


where  I  is  the  NxN  identity  matrix.  In  a  similar  fashion,  the  projection 

i 

operator  on  the  subspace  M  along  the  subspace  M  x^j-^  m]  is 

specified  by 


P  %[i,m]  =  A  %[i,m]  ^A  %[i,m]  A  %[i,m]^  A  %[i,m] 

(6.2.11a) 

It  is  to  be  noted  that  the  following  projection  operator  identity 
holds  as  is  apparent  from  expressions  (6.2.7)  and  (6.2.11a). 


P  %[!,„]  ‘  <6-2'1 

The  complement  of  the  projection  operator  (6.2.11b)  is  formally  given 
by 


^N[i,m] 


=  I  -  P 


^N[i,m] 


(6.2.12) 


A  particular  estimate  x^j-^  of  the  vector  x^  can  be  specified  as 
the  projection  of  x^  on  the  subspace  M  x^^  m]>  c^at  is 


^[i.m]  *  P  %[i,m]  % 


(6.2.13) 


The  error  vector  relative  to  estimate  x^^  and  x^  is  then  given  by 


%[i,m]  "  ^  *  *N[i.m] 


-  Pl 


%[i,m]  % 


(6.2.14a) 

(6.2.14b) 


which  is  expressed  as  a  projection  of  the  vector  x*j  on  the  complement 
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subspace  of  M  ^  j .  It  can  be  straightforwardly  shown  that 

%[i,m]  1  M%[i,m]  f6-2'15) 


where  1  denotes  orthogonality,  that  is,  the  error  vector  is 
orthogonal  to  the  subspace  M  mj  •  The  vector  space  formulation 
described  in  this  section  is  suggestively  depicted  in  Fig.  6.2.1. 


6.3  Linear  Prediction  and  Projection  Operator 

In  this  section  we  will  define  three  methods  of  linear  predictions, 
namely,  forward  prediction,  backward  prediction,  and  delayed  backward 
prediction.  These  methods  will  play  a  central  role  in  the  algorithmic 
solution  procedure  to  be  developed. 


6.3.1  Forward  Prediction 

The  m-th  order  forward  prediction  is  referred  to  as  that 
specific  procedure  for  estimating  the  column  vector  x^  and  by 
means  of  a  linear  combination  of  the  set  of  m  shifted  vectors 
{s'Lxn,  S2^,  ...  ,  3%}  and  {S1^,  S2^,  ...  ,  S™^),  respectively. 
Considering  the  projection  operator  defined  in  Section  6.2,  the 

A  A 

associated  estimates  m]  anc*  m]  are  seen  to  have  the  form 


%[l,m]  *  P  %[l,m]  % 
•^N[l,m]  “  P  %[l,m] 


The  difference  between  the  estimate  x^  and  the  given  vector  x^  is 
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called  the  forward  prediction  error  vector  and  is  specified  by 


■%,m  -^N  ^N[l,m] 


(6.3.2a) 


while  the  forward  prediction  error  vector  of  is  of  the  form 


%,m  =  %  %[l,m] 


(6.3.2b) 


Now  these  error  vectors  are  each  orthogonal  to  the  subspaces  M  ^ 

and  M  x^^  m]>  respectively.  Use  of  complement  projection  operators 
defined  by  (6.2.10)  and  (6.2.12)  yields 


-%,m  P  ^N[l,m] 


(6.3.3a) 


%,m  =  P\[l,m]  *N 


(6.3.3b) 


6.3.2  Backward  Prediction 

The  m-th  order  backward  prediction  is  that  procedure  of 
estimating  the  column  vector  Smx^  and  Sm^  by  a  linear  combination 
of  the  set  of  m  shifted  vectors  S^x^,  •••  >  Sm  ^x^}  and 

...  ,  Sm  respectively.  In  the  same  manner  as  with 

forward  prediction,  by  applying  the  projection  operator,  it  can  be 
shown  that  the  backward  estimate  is  given  by 


%[0,m-l]  =  P  %[0,m-l]  S  % 


(6.3.4a) 


where  the  double  caret  notation  designates  backward  prediction.  The 
backward  prediction  error  vector  is  then  found  to  be 
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A  c  m 

■^N,m  P  ^N[0,m-l]  ^  -^N 


(6.3.4b) 


6.3.3  Delayed  Backward  Prediction 

The  m-th  order  delayed  backward  prediction  is  similarly  defined 

to  be  that  procedure  of  estimating  the  column  vector  and 

1  2 

S  ^  by  a  linear  combination  of  the  sets  of  vectors  {S'S^*  S  x^. 


,  and  (S1^,  S2^,  *  Su\^},  respectively.  It  can  be 


shown  that  the  delayed  backward  prediction  is  given  by 

*  r.  ,.111+1 

%[l,m]  ?  %[l,m]  S 

while  the  delayed  backward  prediction  error  is  specified  by 


,x  _  _c  cm+l 

P  %[l,m]  S  % 


(6.3.5a) 


(6.3.5b) 


A  little  thought  will  convince  oneself  that  the  projection  operators 
P  can  be  expressed  as 


P  %[l,m]  =  A  %[l,m][A  ^N[l,m]  A  ^[l.m]3"1  A '^[l,m] 


”0  0.  .  .  .  o” 

■» 

A  %-l[0,m-l]A  %-l[0,m-l] 

A  %-l[0,m-l] 

.  4 

- 

-r' 


0  0  ....  o 


A  %-l[0,m-l] 


(6.3.6) 


J 
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This  formula  is  straightforwardly  obtained  by  making  use  of  the 
structure  of  matrices  A  x^q  and  A  defined  by  (6.2.8)  and 

(6.2.9).  The  relationship  between  the  backward  error  and  the  delayed 
backward  error  is  then  readily  found  to  be 


m 


[0 


£-1..3 


[0 


— N-l,m 


] 


(6.3.7a) 

(6.3.7b) 


It  then  follows  that  the  N-th  delayed  error  is  equal  to  the  (N-l)-st 
backward  error,  that  is 


m  (N)  =  m  (N-l)  (6.3.8a) 

d£  (N)  =  b£  (N-l)  (6.3.8b) 

The  relationship  between  forward,  backward,  and  delayed  backward  is 
suggestively  depicted  in  Fig.  6.3.1. 


6.4  Decomposition  of  Projection  Operators 

The  development  of  a  computationally  efficient  algorithm  is  depen¬ 
dent  on  the  decomposition  of  the  above  projection  operators.  This 
decomposition  makes  use  of  the  specific  matrix  structure  referred  to 
as  shift  invariancy.  A  matrix  which  has  a  displacement  rank  3  will 
possess  this  shift  invariancy  (see  Chapter  5) .  In  this  section,  the 
shift  invariant  structure  is  utilized  to  decompose  projection 
operators.  The  formulae  obtained  in  this  section  will  be  used  for 


Forward  Backward  and  Delayed  Backward  Relationship 
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the  development  of  order  update  recursions  in  section  6.5. 

First,  we  will  discuss  the  decomposition  of  the  projection 
operator  P  x^q  mj*  This  projection  operator  P  x^j-q  m]  raay  be 


expressed  as 


P  %[0,m]  =  A  -\[0,m]  [R  %[0,m]rl  A\[0,m]  (6*4'1) 

which  is  obtained  by  substituting  i  =  0  in  expression  (6.2.7).  The 


matrix  R  x^j-q  is  defined  by 


R  %[0,m]  ”  A  %[0,m]  A  %[0,m] 


(6.4.2) 


Substitution  of  expressions  (6.2.8)  and  (6.2.9)  into  (6.4.2)  yields 


^N[0,m] 


A  ^[l.m]  % 


%  A  ^[l.m] 


(6.4.3a) 


%[l,m] 


where  R  x^^  m"]  is  defined  by  substituting  1  in  place  of  0  in 
expression  (6.4.2).  If  we  denote  the  inverse  of  matrix  R  Xj^  m]  by 


,  it  then  follows  that 


0 . 0 

R  %[0,m]  *  -1 

•  R 


4  A  %[l,m]  R"  %[l.m] 


(6.4.4) 


where  I  denotes  the  m  x  m  identify  matrix.  Upon  examination  of 
expression  (6.4.3a)  and  (6.4.4),  it  can  be  readily  shown  that 

in  R  %[0,m]  =  %  •  ° 

where  _e  denotes  the  (m  +  1)  x  1  unit  basis  vector  and  u^| ^  is  a 
(m  +  1)  x  1  column  vector  given  by* 


%fl 


A  %[l,m]  R  ^j[l,m]^/fN,m 


in  which 


£ 

is  a  scalar  defined  by 

N,m 


,e  +  x 

«  y  £„ 

N,m 


<**.»> 


X 

%,m 


(6 


(6 


In  a  similar  fashion,  let  us  define  a  matrix  R  mj  by 

R  %[0,m]  A  ^%[0,m]  A  %[0,m] 


(6 


It  then  follows  that 

4l  R  %[0,m]  =  i  (6 

where  v  , -  is  a  column  vector  expressed  by 

tffrl 

if!  “  Cl*  “i  A  %[l,m]  R'i[l,m]]/  f5,m  *  (6> 

Taking  the  complex  conjugate  vector  transpose  of  expression  (6.4. 
yields 


.4.5) 


.4.6) 


.4.7) 

.4.8) 

.4.9) 

4.9a) 

9), 


*In  general,  e^  represents  the  standard  unit  basis  vectors  whose 
components  are  also  zero  except  for  its  k-th  which  is  one. 
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R  %[0,m]  Vl  "  % 


(6.4.10) 


The  inverse  of  the  matrix  R 


^N[0,m3 


is  found  to  be 


R  %[0,m] 


,-l 


+  v  .  , ,  u  , . 

N,m  - nH"l  - nH-1 


(6.4.11) 


Substitution  of  expression  (6.4.11)  into  (6.4.1)  then  leads  to  the 
following  relationship. 


P  ^N[0,m]  P  %[l,m]  +  A  ^N[0,m]  ^N,m  -^nri-1  ^nrt-1  A  %[0,m] 


(6.4.12) 


After  a  simple  algebraic  manipulation,  the  projection  operator 
P  Xjjjjo  m]  is  decomposed  by  the  following  relationships 


P  ^N[0,m]  P  ^N[l,m]  +  P  %,m 

=  P  ^[l.m]  +  (I  ~  P  ^NCl.rn]5  P  %,m 


(6.4.13a) 

(6.4.13b) 


+  P  £  -  Cl  -  P 


%[l,m]  +  P  %,m 


%[l,m] 


) 


(6.4.13c) 


where  it  is  readily  shown  that  P  ^  is  a  projection  operator  onto  the 


V 

subspace  spanned  by  along  the  subspace  which  is  orthogonal  to  the 


subspace  spanned  by  m  and  is  defined  by 


P  e 


1  x 


(Si  J 


,m  e  ,m  ,m 
N,m 


(6.4.14) 


^  ‘  ^  - 1  -  *‘‘“t 


Furthermore,  expression  (6.4.13c)  leads  to  the  following  relationship 


1  ~  P  %[0,mfl]  "  P  “  P  %[l,m]^  (6.4.15) 

The  projection  operator  decomposition  as  expressed  in  (6.4.15)  will 
be  used  to  find  a  backward  error  recursion  in  the  next  section. 

Next,  we  will  decompose  the  projection  operator  P  x^j-^  which 

is  necessary  to  compute  the  forward  prediction  error.  The  projection 
operator  P  nH.1j  is  given  by 

P  %[l,iiH-l]  =  A  ^NCl.m+l]  %[l,mfl]^  A  %[l,mfl] 

(6.4.16) 

which  is  obtained  by  substituting  i  *  1  in  expression  (6.2.7).  The 
matrix  R  x^j-^  is  defined  by 


R  %[l,nH-l]  *  A  %[l,m+l]  A  %[l,nrfl] 

Substitution  of  expressions  (6.2.8)  and  (6.2.9)  with  i  *  1  into 
(6.4.17)  yields 


%[l,mfl] 


(6.4.17) 


- 

.7  „nrH 

R  %[l,m] 

A  %[i,«] s  % 

(S’*1  J*>+  A  %[l>m] 

(S^1  s"*1  % 

(6.4.18) 


It  then  follows 
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%[l,nri-l]  R 


,ciiH-l  ,+.  _-l 

0 . 0  (S  %>  A  %[!,«] 


(6.4.19) 

Upon  examination  of  (6.4.18)  and  (6.4.19),  the  following  relationship 
can  be  derived 


-TTT+-1  R  -^N[l,nri-l]  ”  ^nri-1 


(6.4.20) 


where  e  . ,  is  a  unit  basis  vector  whose  nri-ls  element  is  1  and  s  , , 

- tn+1  — m+1 

is  a  (m  +  1)  x  1  column  vector  defined  by 

+  r  t  nri-1  v  +  .  „-l  ,1/  cd 

Vl  "  [-(s  V  A%[l,m]  R  %[!,«]'  1]/  fN,m 

(6.4.21) 


in  which  f„  is  a  scalar  defined  by 
N,m  J 

cd  ,  nrf-1  x+  ,x  ,,y  .  f  ,x 

fN,m  *  (S  V  4,m  =  (%,m> 


(6.4.22a) 


(,y  v  hx  fr 

^N-l,m'  %-l,m  "  Tl-l.m 


(6.4.22b) 


Relationship  (6.4.22b)  is  obtained  from  (6.3.7a)  and  (6.3.7b).  After 
applying  a  similar  analysis  to  the  matrix  R  nrt-l]  can  s*lown 


t:.,  r 


—nrt-l  R  %[l,  nrt-l]  *  ^trrt-1 


(6.4.23) 


where  t  ^  Is  a  (m  +  1)  x  1  column  vector  expressed  as 

— mfl  *  ^  A  R  1^/  fN,m* 


(6.4.24) 


Applying  the  vector  transpose  operation  to  both  sides  of  expression 
(6.4.23) ,  we  have 


■^Nfl.mfl]  — mfl  -—mfl 


(6.4.25) 


The  inverse  of  the  matrix  R  x^^  nrfl]  rea^^y  found  to  be 


+  flf  t  , ,  si, 

N  ,m  —nrfl  —mfl 


(6.4.26) 


Substitution  of  expression  (6.4.26)  into  (6.4.16)  then  yields 


P  %[l,mfl]  =  P  %[l,m]  +  A  %[l,mfl]  ^.m-^irfl  -^mfl 


*  A  %[l,mfl] 


(6.4.27) 


After  a  simple  algebraic  manipulation,  equation  (6.4.27)  is  compactly 
expressed  as 


P  %[l,nrfl]  ’  P  2»[l,m]  +  P  4,m 


(6.4.28a) 


+  (I  * p  p  4,»  (6-4-28M 
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•  P  +  P  <£,«  (I  -  P  %[!,«]>  (6-4-28c> 
x 

where  it  is  readily  shown  that  P  is  a  projection  operator  onto  the 

9  in 

x 

subspace  spanned  by  d^  m  along  the  subspace  which  is  orthogonal  to  the 

y 

subspace  spanned  by  and  is  defined  by 


TJ.m 

Furthermore,  equation  (6.4.28c)  can  be  expressed  in  the  form 

1  '  P  %Cl,n*l]  ‘  (I  '  P  !$,„>  «  -  P  %[!,.])  <6-4-W 

Expression  (6.4.30)  will  be  used  to  find  the  forward  error  recursion 
in  the  next  section. 

In  a  similar  manner,  the  following  relationship  may  be  also 
obtained 


1  “  P  -s[0,nri-l] 


-  (I  -  P  «  -  P  <6.4.31) 


1  '  P  '  (I  -  P  4,»>  «  -  P  %[!.»]>  (6'4-32) 


where  the  projection  operators  P  e"|  m  and  P  ^  are  defined  by 


p  _ _ I —  ey  <ex  y 

fe  *  %,m  v%,m' 
£N,m 


P  4,»  <4,/ 

N,m 


(6.4.33) 


(6.4.34) 


Expressions  (6.4.31)  and  (6.4.32)  will  be  used  to  find  the  recursion 

Jj..- 


y  y 

of  forward  error  zf.  and  backward  error  b„  . 
— w  ,m 


6.5  Order  Update  Recursions 


In  this  section,  we  describe  the  order  update  recursive  formulas 

St 

which  recursively  compute  the  optimum  nrfl  order  prediction  error 
from  the  optimum  m-th  order  prediction  error.  Expressions  (6.4.15), 
(6.4.30),  and  (6.4.31)  and  (6.4.32)  play  a  central  role  in  obtaining 
these  order  update  recursions. 

'Let  us  first  derive  the  order  update  recursion  for  the  forward 
prediction  error  vectors.  Applying  the  projection  operator  (6.4.30) 
to  the  column  vector  yields 


X  f—  —  jX  \  X 

%,m+l  =  (I  "  P  ^I.m^  %,m 


(6.5.1) 


Substitution  of  expression  (6.4.29)  into  this  relationship  then  yields 


%,mfl  “  %,m  "  d  ^N,m  (^N,m)  %,m 


(6.5.2) 


The  order  update  recursion  for  the  N-th  component  of  the  forward 
prediction  error  vector  is  found  to  be 

■  4,»(«  -  7s*2-  <6-5-» 

N-l,m 

where  the  partial-correlation  coefficients  are  specified  by 


r.nri-1 


*N,m  "  %,m  *  <S  V  «  -  P  %[l,m]>  %  (6-5'4) 


In  a  similar  manner,  applying  the  projection  operator  (6.4.32)  to  the 


column  vector 


leads  to 


1 


siUi<N)  ■  y„(«  - 


tc!  ,  (H-l) 

Vl.a 


(6.5.5) 


where 


V»  ■  ^N,m  '  <S"VT  (I  -  P  %[!,»]>  %  <6-5'6) 


Next,  we  will  find  the  order  update  recursion  for  the  backward 


prediction  error  vector.  Applying  the  projection  operator  (6.4.15) 
to  the  column  vector  S^^x^  is  found  to  yield 


,  x  ,T  _  x  .  ,x 

^N.mfl  (I  “  P  4l,r 


(6.5.7) 


Substitution  of  expression  (6.4.14)  into  this  relationship  yields 


^N,m+-1(N)  =  ^N-l,m(N_1)  "  ,e’  -Sn,ih(N) 


(6.5.8) 


where  the  partial  correlation  coefficient  t„  is  specified  by 

N,m 


Sj.rn  =  ^N,m  =  %  (I  ’  P  ^[l,™]^  S  %  (6.5.9) 


Similarly,  applying  projection  operator  (6.4.31)  to  the  column 


vector  S  is  found  to  yield 


since 


u*7  ..  .  1  o__  - 

— N  ,nrf  1  — N-l,m 


(N-l)  - 


-7^  d  _(n) 

fe  * 

N,m 


(6.5.10) 


%m  ~  (4,m)X  4,m  “  %  (I  "  P  %[l,m]>  (6-5-U) 
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£  V 

Next,  we  will  derive  the  recursion  for  f  and  f ,  .  Manipula- 

N,m  N,m 

tion  of  expressions  (6.4.7),  (6.4.28c)  and  (6.4.29)  eventually  leads 
to  the  form 


^.m  CN,m 

rN,mfl  =  N,m  “  r 

N-l,m 


(6.5.12) 


Expressions  (6.4.22b),  (6.3.4b)  and  (6.4.13c)  yield  the  recursive 


formula 


fr  =  fr  SN,m  Sl.m 

rN,rtri-l  N-l,m  "  fe 

N,m 


(6.5.13) 


Consequently,  expressions  (6.5.2),  (6.5.5),  (6.5.8),  (6.5.10), 
(6.5.12)  and  (6.5.13)  represent  the  order  update  recursions. 


6.6  Time  Update  Recursions 


As  a  new  element  of  the  time  series  is  observed,  the  partial 
reflection  coefficients,  forward  errors,  and  backward  errors  may  be 
recursively  computed  by  making  use  of  these  values  obtained  at  the 
last  time  instant.  This  being  the  case,  these  parameters  are  said  to 
be  "time  updated"  for  each  new  data  point. 

The  matrix  A  x^^  may  he  expressed  in  the  recursive  form 


A  %-l[i,m] 


0  0  •  •  •  0 


A  %[i,m]  ”  PN  A  %[i,m] 


(6.6.1) 


where  is  the  N  x  N  projection  matrix  given  by 


101 


P  =  e  e 
N  — N  % 


(6.6  2) 


in  which  e„  is  an  N  x  1  unit  base  vector.  The  matrix  R 
— N 

also  be  expressed  as 


^[i,mj  may 


R  %-l[i,m]  R  %[i,m]  A  -%[i,m]  PN  A  ^N[i,m]  (6.6.3) 


.-1 


It  then  follows  that  the  matrix  R  m]  recursively  updated 

by  (see  Appendix  C) 


-1 


,-l 


.-1 


R  %-l[i,m]  R  %[i,m]  +  1  -  y.  M  R  %[i,m] 


A  v.. 


■^N[i,m]  PN  A  %[i,m]  R  %[i,m] 


,-l 


(6.6.4) 


where  y.  „  is  a  scalar  defined  bv 
i,m,N 


Yi,m,N  “  %  A  -%[i,m]  R  %  A  %[i,m]  % 


(6.6.5) 


Premultiplying  expression  (6.6.4)  by  (I  -  P„)  A  x  r.  -»  and  then 

N 

4* 

postmultiplying  that  result  by  A'^^,  ^  (I  -  P^)  leads  to  the 
recursive  relationship 


^N-lti.m] 


=  (I  -  P„)  P  > 


P  %[i,m]  (I  “  PN^ 


+  f-Y.  nN  Cl  -  V  P  %[!,»]  PS  P  '  PS] 

1  5  lU  y  jN 


(6.6.6) 
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Since  the  vectors  ^  and  z^.  are  elements  of  vector  space  the 
time  update  recursion  is  given  by 

(I  "  P  ^uCi.ml^  %  "  ^N-l  (I  ~  P  ^N-l[i,m])  %i-l 

“  %  V(1  -  P  4  (6-6*7) 

where  V(I  -  P  x^^  designates  the  time  difference  of  the  projection 
operator  defined  by 

,(I  '  P  %[!,»]>  *  PN  '  P  %[!,»]  + 


(6.6.8) 


0 

P  ^N-l[i,m] 

0 . 0 


Substitution  of  expression  (6.6.6)  into  this  expression  yields 


v{1  -  P  *a[i,«]>  *  PN  -  p  ^[i.m]  +  (1  -  V  p  %[i.»]  (I  -  V 


1  ”  yi,m,N 


-d  -  PN)  P  PN  p  %[i,m]  (I  “  PN) 

(6.6.9) 


Expression  (6.6.9)  is  straightforwardly  carried  out  by  a  simple 
algebraic  manipulation  and  yields  (see  Appendix  C) 


V(I  -  P 


%[i,m] 


) 


l  -  y 


i,m,N 


(I-P^[i,m]>  PK  (I"P%[i,^ 


(6.6.10) 
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Expression  (6.6.10)  is  used  to  find  the  time  update  recursion  formula. 
The  partial  correlation  coefficient  m  is  recursively  calculated  by 


N,m 


SN-l,m  + 


4f-l,m 


b.-:  ,  *(N-i>  c  (n) 


.x 

%,nT 


1  -  Y 


l.m.S 


(6.6.11) 


In  a  similar  manner,  the  partial  correlation  coefficients  t^  ^  is 
recursively  calculated  by 


^.m  “  *11-1, m  + 


%y„ 


hx 

%-l,mv 


*(N)  b*  ,  _(N-1) 


1  -  Y 


l,m,N 


(6.6.12) 


The  time  update  recursion  for  forward  error  is  found  to  be 


z*  *(N)  £*,  (N) 

,e  ,  ,m _ — N,m 

N-l,m  1  -  Y-, 


l,m,N 


(6.6.13) 


The  backward  error  is  also  given  by 

r  r  m*(N)  ^ 

fr  _  fr  .  ~N,m  -N,m 

N,m  N-l,m  1  -  y 


o,m-l,N 


(6.6.14) 


A  recursive  formula  for  auxiliary  parameter  Y^  m  can  be  obtained  by 


using  relationship  (6.4.28c)  to  yield 

x 


Y1,uH-1,N  =  Yl,m,N 


(6.6.15) 


N-l,m 


Finally,  y  can  be  computed  by  using  the  following  relationship 

O )  IQ)  IN 


£  m(N)  *(N) 
^  ,  ~N,m _ — N ,m _ 

Yo,nH-l,N  Yl,m,N  fe 

N  ,m 


(6.6.16) 
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which  is  directly  obtained  from  expression  (6.4.13c). 

Thus  we  can  use  equation  (6.6.11)  and  (6.6.12)  to  update  the 
partial  correlation  coefficients.  Equations  (6.6.13)  and  (6.6.14)  can 
be  used  to  time  update  the  forward  and  backward  covariance  errors 

£  X 

£,  and  flY  .  The  auxiliary  parameters  y,  „  and  y  are 

N,m  N,m  r  l,m,N  o,m,N 

recursively  computed  by  expression  (6.6.15)  and  (6.6.16),  respectively. 


6.7  An  Algorithm  for  Recursive  ARMA  Spectral  Estimation 


In  this  section,  we  summarize  the  recursive  ARMA  spectral  estima¬ 
tion  algorithm  developed  in  the  previous  sections.  For  programming 

X  V 

convenience,  the  following  notations  shall  be  used:  e  (m) ,  e^On), 

n  n 

b*(m),  by(m),  f^(m),  f^(m) ,  s  (m) ,  t  (m) ,  y  (m)  and  y  (m)  in 
n  n  n  n  n  n  o,n  l,n 

place  of  e*  (N) ,  J.  (N) ,  bX  (N) ,  b^  (N) ,  fft  ,  f*  ,  sXT  ,  tM  , 

^  — N,m  -N,m  -N,m  ’  -N,m  N,m  N,m  N,m  N,m 

y  „  and  y  ,  respectively.  At  each  new  data  point,  the  para- 
0, m,N  l,m,N  rt  r 

meters  are  recursively  time  updated  (see  section  6.6)  and  order  updated 
from  m=0  to  m*p-l  (see  section  6.5).  The  recursive  ARMA  spectral 
estimation  algorithm  can  be  presented  as  follows. 


Step  1  Initial  Condition  (Time  Update  n*l) 
e*(0)  =  ey(0)  =  bX(0)  =  by(0)  *  0 

f1(0)  -  f£(0)  *  x^y*,  si(0)  -  ti(0)  *  0  for  i  *  0,  ...  ,  p-1 

Step  2  Initial  Condition  (Order  Update,  m*0) 

eX(0)  =*  bX(0)  -  x  ,  cy(0)  -  by(0)  =■  y 
n  n  n  n  n  n 
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Yl,n-l("1)  '  =  Yi,n(0>  “  0  i  -  0,1 

fn(0)  =  fn(0>  *  fl  i CO)  +  y*  x 
n  n  n-1  n 

Step  3  Order  Update  Recursions 

(m  =  0,  1,  ...  ,  M  for  M  =  min  (p-1,  n-1)) 

(i)  Forvard  Error 

x  „  s  (m) 

=  £n^  -  bn-l(m> 

fn-l(m) 


V/  v  t*(n) 

en(nrfl)  =  £n(m)  “  7T~77*  bn-l(m) 
fn-lCm) 

(ii)  Backward  Error 

v  „  t  (m) 

b  (nri-1)  -  b  <m)  -  -r -  e*(m) 

n  X  f£(m)  n 

n 


y,  v  s*(m) 

bl(nrt-l)  *  by  (m)  -  -2 - 

n  n-1  e .  . *  _n 

f  (m) 
n 

(iii)  f  (m) ,  f^(m)  and  7,  (m) 

11  n  1 ,  n 

e  e  s  (m)  t  (m) 

flCnH-l)  =  f£(m)  -  -2 - a - 

n  n  r 

fn-l<m> 


ey(®) 


if  n  <_  p 


r  r  s  (m)  t  (m) 

f>D  »  <_,(.)  -  -a--  ^  -  if  n  <  p 


f>> 


b*  . (m)  by  .. (m) 

r,  (nrf-1)  -  y  (m)  +  -2^i - Sli - 

■L.n  l.n  r 

*  /  \ 


fn-l(m) 
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£X(m)  e^(m) 

Y.  (nrt-1)  =  Y,  (m)  +  — — - - - - 

0,n  l,n  ,e 

£N,m 


Step  4  Time  Update  Recursions  (m  =  0,  1,  ...  ,  M) 


(i)  Partial  Correlation  Coefficients 


(m)  £X(m) 

s  (m)  »  s  .  (m)  +  — r~ - 7~t - 

n  n-1  1  -  Y,  _(m) 


'l.n' 


e^(m)  *  bX  (m) 

t  (m)  *  t  (m)  +  — t - 

n  n-1  1  -  Y,  _(n) 


'l,nv 


(ii)  fe  (m)  and  fr  , (m) 
n-1  n-1 


f>) 


e^(m)  *  £X(m) 

f  .  (m)  +  —S - — 

n-1  1  -  Y-,  _(m) 


'l.n' 


if  n  >  p 


f>> 


b^(m)  bX(™) 

f  .  (m)  +  —2 - B__ 

n-1  1  -  Y«  _(m) 


0,n 


if  n  >  p 


Step  5  Let  n  *  n+1,  if  n  j*  N  go  to  Step  2 


Step  6  End  of  Algorithm 


In  above  N  is  taken  as  a  time  index  of  a  pair  of  the  last 
observations  and  y^. 
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6.8  Numerical  Examples 

To  test  the  recursive  ARMA  spectral  estimation  algorithm,  the  time 
series  expressed  by  (4.5.2a),  (4.5.2b)  and  (4.5.2c)  were  generated. 

A  program  listing  of  the  fast  algorithm  used  in  obtaining  the 
denominator  coefficients  of  the  ARMA  model  is  illustrated  in  Appendix 
D.2.  As  a  first  example,  64  data  samples  were  generated  according  to 
expressions  (4.5.2a),  (4.5.2b)  and  (4.5.2c).  These  data  samples  are 
plotted  in  Fig.  6.8.1(a).  The  fast  algorithm  was  then  applied  to 
this  64  observations  to  obtain  an  ARMA  spectral  estimate  with  model 
order  (4,4).  The  forward  error  sequence  cx  ,(n)  (n  *  1,  ...  ,  64) 

n ,  4 

is  plotted  in  Fig.  6.8.1(b).  Comparing  Figures  6.8.1(a)  and  6.8.1(b), 

the  forward  error  sequence  is  observed  to  be  more  random  (uncorrelated) 

than  the  given  data  samples  indicating  a  desired  whitening  effect. 

The  resultant  spectral  estimate  is  shown  in  Fig.  6.8.1(c).  The 

resolution  of  the  two  peaks  is  evident,  however,  the  estimated  level 

of  the  first  peak  is  lower  than  that  of  second  peak.  Next,  500  data 

samples  of  the  same  time  series  expressed  were  generated.  These 

samples  are  plotted  in  Fig.  6.8.2(a).  The  forward  error  sequence 

ex  ,(n)  (n  ■  1,  ...  ,  500)  obtained  by  the  fast  algorithm  is  plotted 

~ n,  4 

in  Fig.  6.8.2(b).  It  is  observed  that  the  forward  error  sequence  con¬ 
verges  in  a  relatively  rapid  manner.  In  Fig.  6.8.2(c),  the  resultant 
spectral  estimate  of  model  order  (4,4)  is  illustrated.  The  resolution 
of  the  two  peaks  is  again  evident.  In  addition,  the  height  of 
the  two  peaks  are  equal  as  desired.  As  these  examples  illustrate, 
the  fast  algorithm  maintains  a  high  quality  of  spectral  performance. 
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6 . 9  Summary 

A  recursive  algorithm  has  been  proved  and  stated  for  efficiently 
updating  partial  reflection  coefficients  of  an  ARMA  spectral  estimation 
model.  The  computational  requirement  for  the  order  update  recursions 
and  time  update  recursions  are  12(M  +  1)  and  10(M  +  1),  respectively, 
where  M  is  taken  to  be  the  minimum  of  either  p-1  or  n-1.  Numerical 
examples  show  that  implementation  of  premodification  will  result  in 
only  a  small  degradation  of  spectral  performance.  If  q=*0,  the  ARMA 
model  is  converted  to  the  AR  model.  A  recursive  AR  algorithm  can  be 
also  developed  based  on  a  less  general  vector  space  approach  discussed 
in  this  chapter  (see  Lee  and  Morf,  1980). 


.  ■*.  ....  . -  V'T"~ 


Chapter  7 


CONCLUSION 


The  development  of  computationally  fast  algorithms  for  high 

performance  AKMA  spectral  estimation  was  presented.  The  required 

2 

computation  for  the  unmodified  method  was  reduced  to  0(4p  )  by  using 
a  generalized  Levinson's  approach.  Methods  of  data  modifications  were 
applied  to  reduce  the  computational  complexity.  Modifications, 
referred  to  as  post-modification  with  p  =  q  and  pre-  and  post-modifi¬ 
cation,  achieved  a  computational  complexity  of  0(p  log  p) .  A  fast 
recursive  algorithm  with  a  computational  complexity  of  0(p)  was 
developed  based  on  the  pre-modification  method. 

The  spectral  performance  of  these  methods  was  compared  for 
various  numerical  examples.  Spectral  degradation  had  been  expected, 
because  of  the  restriction  t  =  p  an  :  the  underlying  data  modification, 
however,  these  numerical  examples  illustrated  only  a  small  degradation 
in  spectral  performance.  Moreover,  the  spectral  estimation  performance 
of  these  new  methods  has  been  found  to  be  typically  far  superior  to 
such  contemporary  approaches  as  the  Box-Jenkin^  and  maximum  entropy 
methods. 

Finally,  considering  the  above  two  aspects,  namely,  fast  computa¬ 
tional  implementations  and  high  performance  spectral  estimations, 
these  new  methods  promise  to  be  primary  spectral  estimation  tools. 
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Appendix  A 


RECURSIVE  AR  SPECTRAL  ESTIMATION 


A.  1  Introduction 

In  many  relevant  signal  processing  applications,  one  seeks  to 
characterize  the  spectral  density  of  a  time  series  based  upon  a  finite 
set  of  time  series  observations.  Without  loss  of  generality,  this 
sample  observation  set  is  taken  to  be  the  contiguous  set  of  N  real 
valued  measurements  as  given  by 

x(l),  x(2) ,  ...  ,  x(N)  (A. 1.1) 

One  of  the  most  widely  used  spectral  estimation  models  is  obtained  by 
postulating  the  following  autoregressive  (AR)  structure 

x(n)  +  a.x(n-l)  +  —  +■  a  x(n-m)  =  e(n)  (A. 1.2) 

i  m 

in  which  s(n)  is  a  white  noise  time  series  with  zero  mean  and  variance 
2 

0£  .  Our  object  will  be  that  of  modeling  an  underlying  time  series 
(x(n)}  with  the  AR  model  structure  (A. 1.2)  in  which  the  coefficients 
are  estimated  from  the  given  finite  set  of  observations  (A. 1.1).  This 
is  readily  achieved  by  applying  the  well  known  one-step  predictor. 

An  m-th  order  one-step  predictor,  by  definition,  estimates  the 
value  of  a  random  time  series  using  a  linear  combination  of  the  most 
recent  m  samples.  Namely,  the  sample  x(n)  is  estimated  by  means  of 

0 

the  relationship 


122 


m 

x(n)  »-Ia,  x(n-k) 
k-1  * 


(A. 1.3) 


The  difference  between  this  predicted  value  and  the  observed  value 
x(n)  over  the  observation  interval  is  called  the  prediction  error 
and  is  specified  by 


e(n)  -  x(n)  -  x(n) 


m  <  n  <  N 


(A. 1.4) 


or 


m 


e(n)  =  x(n)  +  £  a,  x(n-k) 

k=l  k 


m  <  n  <  N 


(A. 1.5) 


Writing  these  error  expressions  in  matrix  form  yields 
_e  ■  x  +  Xa 

where  a,  e and  x  are  m  x  1,  (N-m)  x  1,  and  (N-m)  x  1  column  vectors, 
respectively,  given  by 


a  “  [ar  ...  ,  a J1 


e  ”  [e(nri-l),  e(nr+-2),  ...  ,  e(N)] 


x  *>  [x(nrt-l) ,  x(nH-2) . x(N)] 


(A. 1.7a) 
(A. 1.7b) 
(A. 1.7c) 


4 
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where  the  superscript  T  denotes  the  transpose  operation. 

The  a^  coefficients  are  to  be  now  selected  so  as  to  cause  each 
of  the  predictor  error  terms  e(n)  to  be  close  to  zero.  This  selection 
process  will  give  rise  to  the  so-called  optimal  one-step  predictor. 

To  achieve  the  required  objective  of  setting  the  e(n)  to  be  near  zero, 
one  typically  appeals  to  the  least  squares  method  which  minimizes  a 
squared  error  criterion  of  the  form 

f(a)  =  eT  We  (A. 1.8) 

where  W  is  an  (N-m)  x  (N-m)  nonnegative  definite  square  matrix.  The 
minimization  of  this  quadratic  functional  with  respect  to  the  column 
vector  a  is  straightforwardly  carried  out  and  results  in 

XT  WXa°  =  XT  Wx  (A. 1.9) 


It  can  be  shown  that  the  resulting  power  spectral  density 


estimate  of  the  time  series  {x(n)}  is  then  given  by  (Haykin,  1979) 

2 

o 

(A. 1.10) 


Sx(w) 


[  1  +  a°  e~ja)+a°  e"2ju  +  ...  +  a° 


where  the  a£  coefficients  are  obtained  upon  solving  relationship 

(A. 1.9).  Generally  the  solution  of  relationship  (A. 1.9)  requires  on 
3  3 

the  order  of  m  (i.e.0(m  ))  number  of  multiplications  and  additions 
if  that  relationship  is  directly  used.  This  computational  requirement 
can  be  excessive  in  many  real  time  applications.  It  has  been  recently 
shown  by  Lee  and  Morf  (1980)  that  this  computational  requirement  can 
be  reduced  to  0(m)  by  slightly  reformulating  the  matrix  X  and  column 


i  i  ■la'ini  n  V  in 


vector  x.  In  many  interesting  cases,  fortunately,  the  solution  to 
this  modified  system  of  equations  will  be  close  to  that  of  the  desired 
solution  as  represented  by  expression  (A. 1.9).  In  this  Appendix  the 
method  which  is  identical  to  the  LMS  algorithm  of  Lee  and  Morf  (1980) 
is  presented  with  more  emphasis  on  insightful  development. 

This  general  modification  methodology  shall  herein  be  referred 
to  as  data  modification.  Applying  the  specific  data  modification 
method  referred  to  as  prewindowing,  the  matrix  X  is  reformulated  as  the 
N  x  m  matrix  given  by 

r  •»  T 

0  x(l)  x(2)  .  .  .  x(m)  .  .  .  x(N-l) 

0  0  x(l)  .  .  .  x(m-l)  .  .  .  x(N-2) 

X  -  .  .  0 

•  • 

•  •  #  •  •  • 

•  •  •  •  •  • 

0  0  0...0  x( 1)  .  .  .  x(N-m)  (A. 1.11) 

•  J 

while  the  N  x  1  column  vector  x  is  specified  by 

x*  [x(l),  x(2),  ...  ,  x(N)]T  (A. 1.12) 

If  these  new  entrants  are  substituted  into  relationship  (A. 1.9),  an 
efficient  solution  procedure  for  a0  is  possible.  The  structure  of  this 
reformulated  matrix  X  and  the  column  vector  x  enables  us  to  obtain  a 
recursive  least  square  spectral  estimation  algorithm  which  has  an 
excellent  convergence  behavior  and  a  fast  parameter  tracking 
capability  relative  to  the  former  structure.  The  development  of 
this  algorithm  is  predicated  on  the  utilization  of  projection  operator 
theory  (Luenberger,  1969) .  In  the  sections  which  follow  the  necessary 
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projection  operator  theory  to  be  used  in  the  algorithm  is  described. 

A. 2  Vector  Space  Formulation 

In  this  section,  the  given  spectral  estimation  problem  will  be 
cast  into  a  convenient  vector  space  setting.  It  will  be  assumed  that 
the  following  observations  of  the  time  series  {x(n)>  as  specified  by 

x(l),  x(2) ,  ...  ,  x(N)  (A. 2.1) 

are  given.  This  in  turn  will  give  rise  to  the  associated  column  data 
vector 

Xjj  -  [x  ( 1) ,  x(2)  ,  ...  ,  x(N)]T  (A. 2. 2) 

The  vector  3^  lies  in  the  product  space 

^  -  Rx  R:<  . . .  xR  =  RN  (A. 2. 3) 

This  vector  space  can  be  made  into  an  inner  product  space  by 
defining  the  following  inner  product  between  any  two  elements  x^, 

%  £  ^ 

T  N 

"  V  "  "  %  %  "  x(n) 

n=l 

The  corresponding  induced  norm  of  x^j  is 

ii  \ 

u  n*l 


y(n) 


(A. 2. 4) 


then  given  by 


x  (n) 


"1  *5 


(A. 2. 5) 


i . 
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We  next  define  the  shift  matrix  S  which  is  represented  by  the  N  x  N 
matrix 


(A. 2.6) 


Applying  the  shift  matrix  m  times  to  the  column  vector  x^  is  seen  to 
yield 

s\  3  [0,  ...  ,  0,  x(l),  ...  ,  x(N-m-l),  x(N-m)]T  (A. 2. 7) 
m  zeroes 


We  next  construct  the  subspace  Mx^j-^  which  is  spanned  by  the  set  of 

gmY 


vectors  S 
by 

^[i.m]  *  {S%*  si+% . S\} 


This  subspace  will  be  suggestively  denoted 


(A. 2. 8) 


where  the  first  integer  index  i  may  take  on  any  value  in  the  set 

{0,  1,  ...  ,  m>.  Next,  we  let  Px^^  m]  designate  the  projection 

i 

operator  onto  the  subspace  along  the  subspace  Mx^j-^^j.  This 

projection  operator  can  be  shown  to  have  the  form 


^^[i.m]  *  ^%[i,m]  j^N  i,m  ^N[i,m]_ 


^[i.m] 


(A. 2. 9) 


where  Xx^^  mj  is  the  N  x  (m-i+1)  matrix  composed  of  the  following 
ordered  set  of  column  vectors 


^.[1,.]  ■  si+1s, . SV 


(A. 2. 10) 
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Similarly,  Che  projection  operator  on  the  orthogonal  complement  of 
subspace  Mx^^  is  denoted  by 


i 

P2n(i,»r 1  ■ 


(A. 2.11) 


where  I  is  the  N  x  N  identity  matrix.  It  then  follows  that 

P%[i,m]  %  *  %  if  %  e  %[i,m]  (A. 2. 12) 

i 

P%[i,m]  4  1  ^  ^  (A. 2. 13) 

Expression  (A. 2. 12)  and  (A. 2. 13)  specify  those  properties  of  the 
projection  operators  which  will  be  utilized  when  developing  a  recur¬ 
sive  least  square  algorithm  in  the  next  section. 


A. 3  Linear  Prediction  and  Projection  Operator 

In  this  section,  we  will  define  three  methods  of  linear 
prediction,  namely,  forward  prediction,  backward  prediction,  and 
delayed  backward  prediction.  These  projection  operators  will  play  a 
central  role  in  the  algorithmic  solution  procedure  to  be  developed. 

A. 3.1  Forward  Prediction 

The  m-th  order  forward  prediction  method  is  referred  to  as  that 

specific  procedure  for  estimating  the  column  vector  x^  by  means  of  a 

1  2 

linear  combination  of  the  set  of  m  shifted  vectors  {S  x^,  S  x^,  ... 
Smx^}.  It  then  follows  that  the  m-th  order  forward  prediction 
estimate  of  x^  is  of  the  form 
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while  the  associated  forward  error  vector  is  specified  by 
•%,m  ^N[l,m] 

•  %  +  “A, 

k*l 

k 

Upon  examination  of  the  structure  of  the  shifted  vector  S  x^(k  =  1, 

...  ,  m) ,  expression  (A. 3. 2b)  leads  to  the  aforementioned  prewindowing 
formula  where  X  and  x.  are  given  by  (A. 1.11)  and  (A. 1.12),  respectively. 

The  problem  at  hand  is  to  then  find  the  scalar  constants  a^, 

...  ,  a^  which  minimize  the  squared  forward  prediction  error 

«£>  •  f!i%  -  aN[1.„]ll 2  <A-3-3> 

According  to  the  projection  theorem  (Luenberger,  1969),  f (a)  is 
minimized  when  the  error  vector  is  orthogonal  to  each  of  the  one¬ 
dimensional  subspaces  spanned  by  S*x^(i  =1,  ...  ,  m) .  Thus,  we  have 
the  orthogonality  relationship  expressed  by 

-  2fl[1>m])  -1-  sl2%  for  i  =  1,  2,  ...  ,  m  (A. 3. 4) 

which  takes  the  inner  product  format 

<  %  -  S±%  >  *  0  for  i  *  1,  2,  ...  ,  m  (A. 3.5) 

Substitution  of  expression  (A. 3.1)  into  (A. 3. 5)  yields  the  set  of 
linear  algebraic  equations 


(A. 3. 2a) 

(A. 3. 2b) 
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*  <SV  S%  =*  \  =  _<V  S%> 

k=l 


(A. 3.6) 


for  1=1,2,  ...  ,  m 

for  the  optimum  set  of  a^  prediction  coefficients.  These  equations 
are  called  the  normal  equations  and  can  be  put  into  the  matrix  form 


where 


T  T 

^[l.m]  X%[l,m]  -  =  _X^N[l,m]% 


^[l,m]  =  s2% . 


(A. 3. 7a) 


(A. 3. 7b) 


^al’  a2’  ’  am^  (A.  3.  7c) 

Solving  equation  (A.  3. 7a)  for  a  and  substituting  this  solution  into 

expression  (A. 3.1)  then  yields  the  optimum  prediction  vector 

-1 
i 

T 


%[l,m]  ^N[l,m] 


^[l.m]  ^[l.m] 


^[l.m]^  (A’3’8) 


Upon  examination  of  the  projection  operator  (A. 2. 9)  and  this  expression, 
^N[l  m]  seetl  to  comPactly  specified  by 


~%[l,m]  P— N[l,m]  ~N 


(A. 3. 9) 


Thus,  we  see  that  x^^  is  obtained  by  projecting  x^  onto  the  sub¬ 
space  Mx^^  and  the  ro-th  order  forward  prediction  error  vector  is 
obtained  by  projecting  x*j  onto  the  orthogonal  complement  of 
in  the  i^,  that  is 


%,m  “  P%[l,m]  % 


(A. 3. 10) 


v 


I 

t 


The  corresponding  minimum  mean  squared  error  is  then  defined  to  be 


^N,m  ^N,m  -%,m  ^N,m 


A. 3.2  Backward  Prediction 

The  m-th  order  backward  prediction  method  is  that  procedure  of 


estimating  the  m-th  shifted  column  vector  Smx^  by  a  linear  combination 
of  the  set  of  shifted  vectors  {S^x^,  S^x^,  »  Sm  ^x^}.  This  back¬ 

ward  estimate  is  then  of  the  form 

m-1 


%[0,m-l]  *  "kf0  bkS  % 


(A. 3. 12) 


and  the  backward  error  vector  is  defined  by 


,m 


-N,m  ^N[0,m-l3 


(A. 3. 13) 


In  the  same  manner  as  with  forward  prediction,  by  applying  the 
projection  theorem  it  can  be  shown  that  the  backward  estimate  is 
given  by 


A  p  cm 

^%[0,m-l]  ^N[0,m-ll 


(A .  3 .  '  - ) 


The  backward  prediction  error  vector  is  then  found  to  be 

1 

Vm  *  P^[0,m-1]  S% 

and  the  corresponding  minimum  mean  squared  error  is  obtained  by 


i  A. 3.15) 


eb  UT  V  uT  cm 

f„  *  b.T  b,r  *  b.,  S  x,, 
N ,  m  — N ,  m  — N ,  m  — N ,  m  — N 


(A. 3.16) 
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A. 3. 3  Delayed  Backward  Prediction 

The  m-th  order  delayed  backward  prediction  method  is  similarly 

m  j  ^ 

defined  to  be  that  procedure  of  estimating  the  column  vector  S  x^ 

1  2 

by  a  linear  combination  of  the  set  of  vectors  (S  S  •••  » 

Smx^}.  It  can  be  shown  that  the  delayed  backward  estimate  is  given  by 

„nrt-l 


%[l,m]  "  P%[l,m]  5  % 


(A. 3. 17) 


and  the  delayed  backward  error  is  obtained  by 

1 

4,m  “  P%[l,m] 

The  corresponding  minimum  mean  squared  error  is  measured  by 


(A. 3. 18) 


cd  .T  .  ,T  _m+l 

^N,m  ^N,m  4l,m  =  4l,m  S  % 


(A.  3. 19) 


A  little  thought  will  convince  oneself  that  the  projection 
operation  can  be  expressed 


as 


P%[l,m]  *  ^[l,m] 


^N[l,m]  ^%[l,m] 


^[l.m] 


00.  .  .o' 

r  t* 

~o_o_L  ..o' 

^N-l[0,m-l] 

£%-l[0,m-l]  ^N-l[0  ,m-l] 

^-lCo.m-l] 

The  relationship  between  the  backward  prediction  error  and  the  delayed 
backward  prediction  error  is  then  readily  found  to  be 
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d„  „  -  [°:  K  l1 

■tJ,!!!  •  — N-l,m 


(A. 3. 21) 


,th 


It  then  follows  that  N  delayed  prediction  error  is  equal  to  the 

St 

(N-l)  backward  prediction  error 


d*.  (N)  -  ^  .  (N-l)  (A. 3.22) 

■^N,m  — N-l,m 

The  relationship  of  forward,  backward,  and  delayed  backward  is 
suggestively  depicted  in  Figure  A.l. 


A. 4  Decomposition  of  Subspaces 

The  development  of  a  computational  efficient  algorithm  is 
dependent  on  the  decomposition  of  subspaces.  Subspaces  may  be  decom¬ 
posed  by  appealing  to  the  well  known  projection  theorem  (Luenberger, 
1969) .  The  formulae  obtained  in  this  section  will  be  used  for  the 
development  of  order  update  recursions  in  Section  A. 5. 

Since  the  forward  prediction  error  lies  in  the  subspace 

-N,m 


^N[0,m]  buC  is  orth°g°"al  c°  *%[!,„]»  we  can  ex?ress  ^[0,m] 

the  direct  sum  of  t  and  (e.,  },  that  is 

-N,ra 


as 


^N,m[0,m]  ^f[l,m]  ®  -%,m 


(A. 4.1) 


where  (e^  denotes  the  subspace  spanned  by  the  forward  prediction 

error  vector  .  The  projection  operator  on  the  subspace  (e„  }  is 

— N,m  r  -N,m 

defined  bv 


P%,m  "  %,m  (%,ra  %,m'>  ^N,m 


(A. 4. 2) 
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Relationship  (A. 4.1)  can  be  readily  shown  to  yield  the  following 

decomposition  of  the  projection  operator 
i  1 


P%[0,m]  (I  ”  P-% ,m^  P%[l,m] 


(A. 4. 3) 


Similarly,  since  the  delayed  backward  prediction  error  d^  ^  lies 
in  the  subspace  m+l]  ^ut  *s  orthogonal  to  Mx^^  m]»  we  obtain 


nri-l]  =  ^[l.m]  ®  x4j,m; 


(A. 4. 4) 


where  {d^  denotes  the  subspace  spanned  by  the  backward  prediction 
error  vector.  The  projection  operator  on  the  subspace  {d^  is 
defined  by 


P— N  ,m 


^N,m  ^N,m  ^N,m^ 


(A.4.5) 


Relationship  (A. 4. 4)  is  found  to  yield  the  following  decomposition  of 
the  projection  operator 


1  1 

P%[l,m4-l]  (I  P— N , m^  P%[l,m] 


(A. 4.6) 


A. 5  Order  Update  Recursions 

In  this  section,  we  describe  the  order  update  recursive  formulae 
which  recursively  compute  the  optimum  nrfl-st  order  prediction  error 
from  the  optimum  m-th  order  prediction  error.  Expressions  (A. 4. 3) 
and  (A. 4. 6)  play  a  central  role  in  obtaining  these  order  update 
recursions. 


V 


135 


Let  us  first  derive  the  order  update  recursion  for  the  forward 
prediction  error  vector.  Applying  the  projection  operator  (A. 4. 6) 
to  the  column  vector  x„  yields 


%,m+l  ^  P^N,m^  %,m 


(A. 5.1) 


Substituting  expression  (A. 4. 5)  into  this  relationship  then  yields 


-  d  /jT  d  )_1  dT 
-^N,nH-l  -^N,m  -^N,m  ^N,m  — N,m  -^N,m 


(A. 5. 2) 


Recalling  expression  (A. 3. 22),  the  order  update  recursion  for  the  N-th 
forward  prediction  error  is  found  to  be 

%,nrfl(N)  *  "  ^N.mfl^H-l.m5  ^N-l,m(N_1)  (A. 5. 3) 

where  the  partial-correlation  coefficients  are  specified  by 

1 

VnH-1  "  4i,N-N,m  *  %  P^j[l,m]  S  %  (A.5.4) 


Expression  (A. 5.1)  leads  to 


T 

•^N,nri-1  ^.m+l 


%,m  ^  P^N,m^  %,m 


(A. 5. 5a) 


The  recursion  for  the  forward  minimum  mean  square  error  is  similarly 
found  to  be 


N,nri-1  N,m 


V 


m+1 


(fb  )"1 
vrN-l,m; 


VnH-l 


(A. 5. 5b) 


Expressions  (A. 5. 3)  and  (A. 5. 5b)  constitute  the  order  update  recursion 
formulae  for  the  forward  prediction. 
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Next,  we  will  find  the  order  update  recursion  for  the  backward 
prediction  error  vector.  Applying  the  projection  operator  (A. 4. 3)  to 
the  column  vector  S^^x^  is  found  to  yield 

b.,  =  (I  -  Pc  )  d  (A. 5. 6 

-N,mfl  — N,m  ^N.m 

Substituting  expression  (A. 4. 2)  into  this  relationship  results  in 


-^N,nrt-1  -^N,m  -^N,m  ^N,m  -^N,m  ^N,m  (A. 5. 

The  order  update  recursion  for  the  N-th  backward  prediction  error  is 
then  specified  by 


V«nW) 


^N-l,m(N_1)  ~  ^N,mfl  (fN,m^  %,m(N) 


(A.  5. 8) 


Expression  (A. 5. 6)  leads  to 


^N,nH-l  -^N  ,m+l  ^N,m  ^  ^-%,m^  ^N, 


m 


(A. 5. 9) 


The  recursion  for  f„  is  next  found  to  be 
N,m 


rb  ~ b  .  j  v  *1 

N,nH-l  =  rN-l,m  ~  ^N.nrt-l  UN,m;  ^N.nrt-l 


(A. 5. 10) 


The  order  update  recursion  formulae  for  the  backward  prediction  are 
represented  by  relationships  (A. 5. 8)  and  (A. 5. 10). 


A. 6  Time  Update  Recursions 


As  a  new  element  of  the  time  series  is  observed,  the  partial- 
correlation  coefficients,  forward  least  square  errors,  and  the  backward 
least  square  errors  can  be  computed  recursively  by  using  the  knowledge 
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of  these  parameters  from  the  last  time  instant.  This  being  the  case, 
these  parameters  are  said  to  be  "time  updated"  for  each  new  data 
point.  These  update  recursions  are  obtained  by  utilizing  a  method 
referred  to  as  projection  operator  decomposition. 

For  the  spectral  estimation  problem  considered  here,  we  decompose 

the  projection  operator  Px^^  into  one  that  projects  on  all  past 

observations  and  another  that  generates  the  correction  due  to  a  new 

observation  x(N).  First,  we  define  the  component  projection  matrix 

P  by 
N  y 


PN  "  eNeN 


(A. 6.1) 


where  e^  is  the  N  x  1  unit  basis  vector  expressed  by 


-  [0,  ...  .  0,  l]1 


(1.6.2) 


Let  us  define  the  column  vectors 


2pN  =  Pn  %  =  ^0.  •••  .0.  x(N)]T 

1  1  1 
2pN  3  PN  %  *  [x(l),  ...  ,  x(N-l) ,  0] 

T 


1  i 


(A. 6. 3) 


(A. 6. 4) 


Note  that  ^  ^  ^  ^  and  similarly  for  x^  Ipf}.  The 

projection  of  3^  on  the  subspace  Mx^^  m]  is  now  decomposed  by 
component  projection  matrix  to  obtain 


(A. 6. 5) 


138 


Multiplication  of  (I-P^)  and  the  matrix  Xx^j-^  mj  yields  the 
called  oblique  matrix 


so- 


%[i,m]  (I  ^[i.m] 


(A. 6. 6) 


whose  last  row  is  the  zero  row  vector.  We  define  the  oblique 
projection  operator  to  be 

Q%[i,m]  =  %[i,a]  [C%[i,m]  C%[ifm]]  "^[i.m]  (A-6'7) 

and  its  associated  orthogonal  complement  by 

1 

%.[!,•]  ‘  1  -  <A-6'8) 

Upon  inspection  of  expression  (A. 6. 7),  we  see  that  the  application  of 
the  oblique  projection  operator  to  the  vector  x^  implicitly  possesses 
the  solution  of  the  prediction  coefficients  at  the  N-lst  stage. 

After  simple  algebraic  manipulation,  relationship  (A. 6. 5)  can  be 
expressed  as 

1 

P%[i,m]  %  3  Q%[i,m]  +  P%[i,m]  PNQ%[i,m]  ^  (A’6‘9) 


The  orthogonal  complement  projection  of  x^  can  be  expressed  as 
1  i 

P%[l,m]  %  =  %  "  %j[i,m]  %  “  P^N[i,m]  PN  Q%[i,m] 


(A. 6 . 10a) 


which  can  be  further  developed  to  the  form 


*%[!.■]  %  "  %  •  %<[i,m]  %  -  N  %[i,m]  % 


1  1 
+  P^N[i,m]  PN  ^%[i,m] 


(A. 6. 10b) 


Considering  the  relationships  (A. 6. 3),  (A. 6. 4)  and  (A. 6. 7),  we  obtain 


P%[i,m]  % 


P%-l[i,m]  ^N-l  '  i  1 

srrrm  % 

_ 

(A. 6. 11) 


in4“l  T 

Premultiplying  [S  x^l  on  both  sides  of  expression  (A. 6. 11)  gives 
the  time  update  recursions  of  the  partial  reflection  coefficients 

i  1 

AnrKL,N  =  AnrKL,N-l  +  ^  P%[l,m]  ?N  %j[l,m]  % 


(A. 6. 12) 


where  i  was  taken  to  be  1.  Furthermore,  operation  of  the  component 
projection  operator  on  both  sides  of  expression  (A. 6. 10a)  yields 
11  1 
PN  P%[i,m]  %  =  PN  Q%[i,m]  %  ~  *1^  P%[i,m]  eN®N  %l[i,m]  % 


PN 


Thus  we  obtain  the  relationship 


tt[i,m]  %  1  “  eN  P%[i,m] 


(A. 6. 13a) 


(A. 6. 13b) 


PN  Q%[i,m]  %  *  T-T.  ~  '  PN  P^[i,m]  % 

X  y  IQ  9  N 


(A. 6. 14) 
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where 


Yi,m,N  “  ®N  P%[i,m] 


(A. 6. 15) 


Directly  substituting  (A. 6. 14)  into  (A. 6. 12),  we  see  that 

1  1 

.  A  ,  [S  ^[l.m^N^  ?%[l.m]% 

TT.nrt-l  "  w-l,nrt-l  ^  1  —  y . 

l,m,N 


(A. 6. 16) 


which  simplifies  to  the  form 


N-l,nH-l  1  -  1T1>mj!I 


(A. 6. 17) 


£  r 

Similarly,  the  time-update  for  f„  and  f„  can  be  obtained  as 

N,m  N,m 


f®  =  f®  ,  + 


N,m  N-l,m  1  -  y 


(A. 6. 18) 


l,m,N 


bl  (N) 

fr  m  fr  ,  ~M,m 


N,m  N-l,m  1  -  y 


(A. 6. 19) 


0,m-l,N 


where 


Y0,m-1,N  "  %[o,m-l] 


(A. 6. 20) 


Thus  we  can  use  equation  (A. 6. 17)  to  update  the  partial  reflection 


coefficients.  Equations  (A. 6. 18)  and  (A. 6. 19)  can  be  used  to  update 
forward  and  backward  prediction  errors,  respectively. 
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A.  8  Summary 

A  recursive  algorithm  has  been  presented  for  efficiently  obtaining 
an  autoregressive  (AR)  spectral  estimate.  To  achieve  a  significant 
computational  improvement,  prewindowing  was  applied,  and  projection 
operators  were  utilized  in  the  vector  space  setting.  Normalizations 
of  the  order  and  time  update  algorithm  yields  more  computational 
advantage  than  the  unnormalized  method.  Interested  reader  may  refer 
to  (Lee  and  Morf,  1980)  and  (Friedlandar,  1980). 


Appendix  B 


ADAPTIVE  SPECTRAL  ESTIMATION 


B.l  Introduction 

In  this  chapter,  we  will  discuss  two  adaptive  techniques, 
namely,  the  Widrow-Hoft  algorithm  (Widrow  and  Hoff,  1960)  and  the 
Iterative  LMS  method.  It  is  well  known  that  the  Widrow-Hoff  algorithm 
is  a  recursive  technique  which  updates  parameters  with  the  arrival 
of  each  new  data  sample.  At  each  recursion,  parameters  are  algo¬ 
rithmically  selected  in  a  least  squares  sense.  As  the  number  of  data 
samples  increases,  the  model's  parameters  "may"  converge  to  the  least 
square  solution  which  is  also  known  as  the  Wiener  solution  (Wiener, 
1949) .  Primary  reason  for  utilizing  the  Widrow-Hoff  algorithm  is 
computational  in  nature.  As  each  new  data  point  is  obtained,  only 
0(p)  computations  are  required  to  update  the  model's  parameters. 

The  Iterative  LMS  method  is  a  technique  which  updates  the  solution 

for  the  linear  system  of  equations  which  approximates  the  Wiener 

equations  (Wiener,  1949).  Although  the  number  of  computations  for 

the  Iterative  LMS  method  to  update  parameters  at  every  new  data  point 
2 

is  0(p  ) ,  the  Iterative  LMS  method  gives  the  exact  solution  to  a  given 
linear  system  of  equations.  To  compare  these  two  techniques,  a  number 
of  examples  are  presented. 
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B.2  Widrow-Hoff  Algorithm 

The  analysis  of  an  adaptive  filter  can  be  developed  by  considering 
the  linear  configuration  shown  in  Fig.  B.2.1.  An  adaptive  filter  is 
composed  of  a  tapped  delay  line,  adjustable  weights  and  summers. 

Delayed  signals  which  are  real  valued  are  weighted  and  summed  to 
form  an  output  signal  d(n)  which  designates  an  estimate  for  the 
desired  signal  d(n) .  At  the  n-th  observation,  a  set  of  delayed 
signals  can  be  formulated  in  a  vector  form  ' 

x  =  [x(n  -  1),  x(n  -  2) .  x(n  -  p)]T  (B.2.1) 

— n 

where  x  is  a  pxl  column  vector.  It  is  also  convenient  to  denote  the 
— n  r 

adjustable  weights  at  the  n-th  iteration  by 

h  =  [h (1),  h  (2),  ...  ,  h  (p)]T  (B. 2.2) 

mi  n  n  n 

where  h^  is  a  pxl  column  vector.  The  estimate  of  the  value  of  d(n) 
based  on  the  vector  (B.2.1)  will  be  taken  to  be  the  linear  combination 

d(n)  *  h^  x 
— n  — n 

P 

=  I  h  (k)  x(n  -  k)  (B.2. 3) 

k*l  n 

The  error  between  the  desired  signal  and  the  estimate  at  the  n-th 
sample  is  given  by 

e(n)  =  d(n)  -  h^  x 
~ n  — n 


(B.2. 4) 
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Fig.  B.2.1  Adaptive  Linear  Configuration 
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The  associated  mean  square  error  is  defined  by 


f(h  )  =  E[e  (n) ] 


Substitution  of  (B.2.4)  into  (B.2.5)  is  found  to  yield 


£<V  •  *dd<°>  - 2  4*  \  +  R* 


(B.2.5) 


(B.2.6) 


where  d. . , (0)  is  the  variance  of  the  desired  signal  d(n),  that  is, 
dd 


*dd(0)  =  E^d  (n)^ 


(B.2.7) 


while  r  and  R  are  the  pxl  cross  correlation  vector  and  the  pxp 
dx  x 

covariance  matrix,  respectively,  defined  by 


rdx  “  ^dx(1)»  W25’  »  *dx(P)] 


(B.2.8a) 


■*,(!)•  »xx(0)>  '  •  *  *xx<P  "  2)l 


(B.2.8b) 


(B.2.8c) 
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and  $  (i)  denotes  the  autocorrelation  sequence  of  the  input  signal 

specified  by 

*xx(i)  =  E^x(tl  +  i)  x(n)^  (B. 2 . 8d) 

It  may  be  observed  from  expression  (B.2.6)  that  the  mean-square  error 
is  precisely  a  second  order  function  of  the  weights and  is  visualized 
as  a  parabolic  function  of  the  weight  variables.  The  adaptive 
process  seeks  the  minimizing  weight  variable  selection  by  using  the 
well-known  method  of  steepest  descent. 

In  seeking  the  minimum  mean-square  error  by  the  method  of 
steepest  descent,  one  first  begins  with  an  initial  guess  of  the  model's 
weight  parameters.  The  next  estimate  is  then  obtained  from  that 
estimate  by  making  a  change  in  the  weight  vector  in  the  direction 
of  the  negative  of  the  gradient  vector.  The  gradient  is  obtained 
by  differentiating  expression  (B.2.6)  to  yield 

7£<V  -  -2  +  2  \  iv, 

If  each  change  in  the  weight  vector  is  made  proportional  to  the 
negative  of  the  gradient,  the  method  of  steepest  descent  leads  to 
the  following  recursive  relationship 


h  . ,  =  h  +  pVf (h  ) 
— n+1  — n  — n 


(B.2.10) 


For  a  sufficiently  small  value  of  u,  the  mean-square  error  at  the 
(n  +  l)-st  step  is  approximately  found  to  be 

f (h  )  »  f(h  )  -  2  y2||  Vf(h  )||2  (B.2. 11a) 


(B. 2 . lib) 


2 

where  ||^f(]ln)||-  is  the  positive  scalar  defined  by 

||vf(h  )||2  =  [vf(h  )]T  [Vf(h  )] 

n  ~~n  — n 

It  may  be  observed  from  eq.  (B.2.11a)  that  the  mean-square  error  is 
reduced  with  each  change  of  the  weight  vector.  For  a  proper  choice 
of  u,  it  has  been  claimed  that  this  algorithm  will  converge  to  an 
optimum  point  regardless  of  the  initial  weights.  (Widrow,  1971) 

The  method  of  steepest  descent  requires  the  determination  of  the 
gradient  vector.  In  practice,  the  true  values  of  these  gradients 
are  seldom  available.  To  overcome  this  difficulty,  the  "LMS  algorithm" 
offers  a  practical  procedure  for  implementing  the  method  of  steepest 
descent.  This  algorithm  uses  gradient  estimates  in  place  of  true 
gradient  values.  These  estimates  may  be  "noisy"  (i.e. ,  contain 
errors)  but  the  effect  of  the  gradient-measurement  errors  is  observed 
to  be  small  in  many  practical  applications. 

A  method  of  measuring  gradients  of  the  mean  square  error  which 
does  not  require  squaring,  averaging  or  differentiating  is  now  given. 
The  mean  square  error  fQi^)  may  be  represented  crudely  by  the 
single  sample  e(n),  the  square  of  the  n-th  error  value.  Then  the 
gradient  vector  is  approximated  by 

7f(h  )  -  Ve2(n)  -  -2e(n)  x  (B.2.12) 

— •  — n 


In  order  to  approximate  the  gradient  vector,  the  present  input-signal 
x  and  its  associated  scalar  error  e(n)  are  used.  Upon  taking  an 
expected  value  on  expression  (B.2.12),  expression  (B.2.9)  can  be 


148 


obtained. 

An  adaptation  cycle  will  proceed  with  the  arrival  of  each  new 
input  vector.  From  eqs.  (B.2.10)  and  (B.2.12),  the  adaptation 
procedure  comprising  the  LMS  algorithm  is  completely  represented  by 
(B.2.13). 


e(n)  ■  d(n)  -  h  x 
— n  — n 


h  , ,  *  h  -  2pe(n)  x 
— n+1  — n  — n 


(B.2.13a) 

(B.2.13b) 


Upon  examination  of  expressions  (B.2.13),  we  can  see  that  the 
computational  requirement  is  0(p) .  In  this  algorithm,  the  selection  of 
P  is  also  an  important  factor.  If  p  is  made  too  small,  convergence  is 
slow.  On  the  other  hand,  if  p  is  selected  to  be  too  large,  the  adaptive 
method  may  not  converge.  In  terms  of  selecting  a  best  p,  the 
interested  reader  may  refer  to  (Widrow,  1971;  Luenberger,  1973; 

Huffman  and  Nolte,  1980). 


B.3  Iterative  LMS  Method 

We  will  now  investigate  the  problem  of  how  to  linearly  filter  an 
observed,  wide-sense  stationary,  discrete-time,  random  time  series 
(x(n)}.  Our  primary  Interest  is  to  best  estimate  the  desired  discrete¬ 
time  random  time  series  (d(n)}  In  the  minimum  mean  square  sense. 

The  problem  is  illustrated  in  Fig.  B.3.1.  Our  objective  is  to  find 
the  transfer  function  H(z)  that  minimizes  the  mean  square  error.  We 
assume  that  the  estimate  of  element  d(n)  is  of  the  form 
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P-1 

d(n)  =  Z  h(k)  x(n  -  k)  (B.3.1) 

k-0 

where  h(k)  are  the  filter  weight  elements.  The  estimate  d(n)  is  then 
seen  to  be  a  linear  combination  of  the  most  recent  p  values  of  the 
observation  signal.  The  mean  square  error  is  found  to  be  a  function 
of  filter  weights  h(k)  and  is  specified  by 

f(h)  =  E[{d(n)  -  d(n)}2]  (B.3.2) 

where  h  is  the  pxl  column  vector  defined  by 

h  -  [h(0),  h(l),  ...  ,  h(p  -  1)]T  (B.3.3) 


Substitution  of  expression  (B.3.1)  into  (B.3.2)  and  taking  the  expected 
value  operation  yields 

f(h)  -  r,(0)  -  2  hT  r.  +  hT  R  h  (B.3.4) 

—  a  —  —ax  —  x 

2 

where  r^O)  *  E[d  Cn)] ,  r^x  is  the  pxl  column  vector  whose  k-th  element 
is  given  by  E[d(n)  x(n  -  k)]  for  k  =  1,  2,  ...  ,  p  and  Rx  is  the  pxp 
matrix  whose  elements  are  given  by  Rx(i,j)  =  E[x(n  -  i)  x(n  -  j)] 

(see  eq.  B2.6) . 

The  optimum  filter  weights  vector  is  readily  determined  by  taking 
the  gradient  of  quadratic  functional  (B.3.4)  with  respect  to  h  and 
setting  this  gradient  equal  to  the  zero  vector.  This  is  found  to 
result  in  the  well-known  Wiener  vector  selection  (Wiener,  1949) . 


(B.3.5) 


I 
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Although  this  approach  is  indeed  attractive  and  typically  results  in 
satisfactory  performance,  it  suffers  one  serious  drawback.  Its 
implementation  requires  apriori  covariance  knowledge  which  is  usually 
lacking  in  many  typical  applications. 

In  order  to  achieve  our  object  without  requiring  any  statistical 
information,  we  introduce  an  estimation  error  criterion  defined  by 


fN(h)  =  Z  [d(k)  -  d(k)]2 

a  k“p 


(B. 3.6) 


It  will  be  beneficial  to  represent  this  error  criterion  in  a  vector 
format.  Let  us  define  the  (N  +  1  -  p)  x  1  estimation  error  vector 


x(p-l)  . 


d(p+l)  x(P+l)  x(p) 


x(N)  x(N-l) 


which  can  be  compactly  expressed 


e  *  d  —  X.  h 


x(l)  h(o) 
x(2)  h(l) 


x(N+l-p)  Ih(p-l) 


(8.3.7) 


(B.3.8) 


Using  these  expressions,  the  square  error  criterion  can  be  represented 


yw  •  <4  -  s  i)1  <4  -  X,,  h) 


(B.3.9) 


Minimization  of  the  functional  (B.3.9)  is  straightforwardly  carried 
out  by  setting  the  gradient  7h  f^(h)  equal  to  zero  and  yields  the 
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following  result. 


£  •  &S  \Tl  «*  4 


— N 


(B.3.10) 


In  general  applications,  the  use  of  this  method  is  not  practical  since 

3 


it  requires  on  the  order  of  p  multiplications  to  invert  the  pxp 

,T 


matrix  [x^  X^].  We  will  next  discuss  a  straightforward  procedure  which 
reduce  this  computational  complexity. 

Upon  examination  of  relationship  (B.3.7)  and  (B.3.8),  we  can  see 
that  when  the  new  data  element  x(N  +  1)  is  provided,  the  equation  error 
can  be  updated  by 


%+l  =  -^N+l  "  \+l  - 


(B.3.11) 


—  « 

d(N+l) 

T 

^N+l 

•  - 

(B.3.12) 


where  x^+^  i-s  the  P*1  column  vector  specified  by 
2Sn+1  =  [x(N+l) ,  x(N) ,  ...  ,  x(N-p+2)]T 


(B.3.13) 


It  is  clear  from  relationship  (B.3.10)  that  we  have  to  invert  the 
matrix 


^+1  ^+1^  *  *N  +  %+l 


(B.3.14) 


The  following  recursive  relationship  may  be  used  to  efficiently  update 
the  required  matrix  inverse 


^+1  *8+1^  1  +  xT+1  ^+1  %H  %f-l 


(B.3.15) 


where 


%+l  3  ^  ^N+l 


(B.3.16) 


After  a  few  simple  manipulations,  the  following  recursion  is  obtained 


„o  .0  „  ^ 

^N+l  “  +  .  T  ■%+! 

1  +  %+1  %+1 


(B. 3.17) 


Recursive  relationships  (B.3.16)  and  (B.3.17)  constitute  a  more 
computationally  efficient  method  than  the  direct  approach  (B.3.10). 

2 

It  can  be  shown  that  the  computational  complexity  is  of  the  order  p  . 


B.4  Numerical  Examples 


In  this  section,  we  shall  demonstrate  the  performance  of  two 
adaptive  methods,  namely,  the  Widrow-Hoff  algorithm  and  the  Iterative 
LMS  method.  This  will  be  accomplished  by  investigating  the  time 
series  whose  elements  are  given  by 


x(n)  **  /20"  sin  (0.1  trn)  +  w(n) 


(B.4.1) 


where  w(n)  is  a  white  Gaussian  noise  with  variance  one.  The  normarized 
Weiner  equation  error  can  be  defined  by 


R  h  -  r, 

1  x  -n  -dx1 


(B.4. 2) 


— dx 


where  is  the  pxp  covariance  matrix  of  the  sequence  {x(n)}  and  r^ 
is  the  pxl  cross-correlation  vector  of  the  sequences  (d(n)}  and  (x(n)}. 
The  above  scalar  value  C(n)  yields  a  normalized  measure  of  how  closely 
the  Wiener  equations  are  being  approximated.  All  graphs  except 
Fig.  B.4.4  provide  the  plot  of  normalized  Wiener  equation  error 
referring  to  expression  (B.4.2)  versus  iteration  number  (i.e.,  the 
number  of  observation  data).  The  desired  signal  d(n)  is  specifically 
chosen  to  be  x(n+l) .  This  yields  a  problem  of  predicting  one  step 
into  the  future.  Unless  specified,  the  covariance  matrix  is  initial¬ 
ized  at  15-th  iteration  number. 

It  can  be  observed  from  Fig.  B.4.1  that  the  normalized  Wiener 

equation  error  of  the  Iterative  LMS  method  converges  to  approximately 

zero  after  2300  iterations,  however,  the  Widrow-Hoff  algorithm 

with  y  *  0.001  fails  to  converge.  In  the  Widrow-Hoff  algorithm, 

the  value  of  u  was  next  selected  to  be  .0001  and  .01  in  Fig.  B.4.2 

and  Fig.  B.4.3,  respectively.  As  we  can  see  on  Fig.  B.4.2,  both  of 

the  adaptive  algorithms  converge  reasonably  close  to  zero.  The 

Iterative  DfS  method  converges  faster  than  the  Widrow-Hoff  algorithm. 

Fig.  B.4.3  illustrates  an  example  which  shows  convergence  behavior 

of  the  Iterative  LMS  method  and  nonconvergence  behavior  of  Widrow- 

Hoff  algorithm.  The  normalized  square  error  ]]h  -  hA||/[|hA[|  where 

n  u  u 

hg  is  the  exact  solution  of  the  matrix  equation  (B.3.5)  are  displayed 
in  Fig.  B.4.4.  The  convergence  behavior  of  the  Iterative  LMS  and 
the  nonconvergence  behavior  of  the  Wiener-Hoff  are  evident. 


for  Two 


Fig.  B.4.3  Normalized  Weiner  Equation  Error  for  Two  Adaptive  Algorithms  (p  “  A)  for 


Fig.  B.4.4  Normalized  Squared  Error  for  Two  Aduptive  Algorithms  (p 
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Fig.  B.4.5  and  Fig.  B.4.6  display  the  method  which  employs 
biased  estimates  for  the  approximation  of  covariance  matrix  elements. 
Fig.  B.4.7  and  Fig.  B.4.8  display  the  method  which  uses  unbiased 
estimates  for  the  approximation  of  covariance  matrix  elements.  Both 
the  biased  and  unbiased  methods  converge  to  zero,  however,  the  biased 
method  starts  with  slightly  large  values  of  normalized  Wiener  equation 
error.  Fig.  B.4.9  and  Fig.  B.4.10  show  the  Iterative  LMS  method 
whose  initial  covariance  matrix  is  the  identity  matrix.  Although  the 
normalized  Weiner  equation  error  at  the  early  stage  of  iteration 
number  are  relatively  large,  this  method  also  converged  to  zero. 

Fig.  B.4.11  and  Fig.  B.4.12  display  the  direct  method.  Upon 
examination  of  Fig.  B.4.5  through  Fig.  B.4.12,  the  direct  method  and 
the  method  of  unbiased  estimate  are  found  to  be  the  best,  since  they 
started  with  a  smaller  normalized  error  and  converged  uniformly  to 
zero. 

Comparing  the  Widrow-Hoff  algorithm  and  the  Iterative  LMS  method 
from  the  convergence  viewpoint,  the  Iterative  LMS  method  is  superior 
to  the  Widrow-Hoff  algorithm. 

B .  5  Summary 

Two  adaptive  techniques  are  compared.  From  a  computational 
viewpoint,  the  Widrow-Hoff  algorithm  is  less  burdensome  than  the 
Iterative  LMS  method.  However,  the  comparison  of  Wiener  equation  errors 
indicated  that  the  solution  from  the  Iterative  LMS  method  satisfies 


Wiener  equations  better  than  that  of  the  Widrow-Hoff  algorithm. 


Iteration  Humber  n 


Normalized  Weiner  Equation  Error  for  Two  Adaptive  Algorithms  fp 


Iteration  Number  n 


Iteration  Number 


KJg.  B.4.9  Normalized  Weiner  Equation  Error  for  Two  adaptive  Algorithms  (p 


Number  n 


Normalized  Weiner  Equation  Error  for  Two  Adaptive  Algorithms  (p  -  10) 


Iteration  Number  u 
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Appendix  C 


DERIVATION  OF  EXPRESSIONS  (6.6.4)  AND  (6.6.10) 


In  this  Appendix,  relationships  (6.6.4)  and  (6.6.10)  which  play 
a  central  role  for  the  time  update  mode  are  derived. 


C.l  Derivation  of  (6.6.4) 


Expression  (6.6.3)  can  be  simplified  to  the  form 


St  .  =  R  -  x  y_ 
n-1  n  — n  -ti 


(C.l. la) 


where  the  (m-i+l)xl  column  vectors  and  ^  are  defined  by 


^n  ”  A  %[i,m]  % 


(C.l. lb) 


=  %  A  %[i,m] 


(C.l.lc) 


It  can  be  seen  that  is  expressed  as  a  sum  of  a  nonsingular  matrix 

and  a  rank  1  matrix.  Expression  (C.l. la)  can  be  also  expressed  as 


l"1.  »  R""5  [l  -  a  b  ]-1  R~ 
n-1  n  -  n 


,-4 


-1  n-% 


(C.l. 2a) 


where  the  (m-i+l)xl  column  vectors  a  and  b  are  defined  by 


a  ■  R  54  x 
—  n  — r 


b  »  y  R 
—  -ti  n 


-k 


(C.l. 2b) 
(C.l. 2c) 


Mil  1 II  I4m  At 
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in  which  R  ^  is  a  NxN  matrix  which  satisfies  the  relationship 
_  2 

R  R  *  R  .  We  will  now  make  use  of  the  following  matrix  inverse 
n  n  n  6 

relationship 


[t  -  a  Vi'1  -  1+  (T  a)  .  b+ 


(C.1.3) 


Substituting  (C.1.2b)  and  (C.1.2c)  into  (C.1.3)  yields 


[I  -  a  b+]-1  -  I  + 


1  -h  f  -h 

- r - i - R  x  *■ 

,  t  -1  n  — n  n 

1  -  y_  R„  x 

*ti  n  — n 


(C. 1.4) 


Expression  (6.6.4)  can  be  obtained  by  substituting  (C.1.4)  into 
(C.1.2a)  along  with  expressions  (C.l.lb)  and  (C.l.lc). 


C.2  Derivation  of  (6.6.10) 


To  simplify  the  complexity  of  notations,  let  us  define  the 


following  compact  notations 


P  -  P 


Q  -  Pv 


%[i,m] 


Y  *  Yi,m,N  =  4  P  % 


It  is  readily  shown  that 


(1  -  y)  Q  *  %  (1  -  y)  %  *  Q  (I  -  P)  Q 


PQ  PQ  -  yPQ 


QP  OP  -  yQP 


(C.2. la) 


(C.2. lb) 


(C.2.1c) 


(C.2. 2a) 


(C.2. 2b) 


(C.2. 2c) 


I 
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Equation  (6.6.9)  may  be  expressed  as 


7(1  -  p)  -  Q  -  p  +  (I  -  Q)  P(I  -  Q) 


+  -- j;—  (I  -  Q)  P  Q  P(I  -  Q) 

Using  relationships  (C.2.2b)  and  (C.2.2c),  we  have 

{Q  -  P  +  (I  -  Q)  P(I  -  Q)>  (1  -  Y) 

-  Q  -  PQ  +  yPQ  -  QP  +  yQP  -  yQPQ 


(C.2.3) 


(C.2.4) 


Substitution  of  (C.2.4)  into  (C.2.3)  yields 

7(1  -  P)  =  (i  -  P)  Qd  -  P> 

Expression  (6.6.10)  can  be  obtained  by  direct  substitution  of 
expressions  (C.2.1a),  (C.2.1b)  and  (C.2.1c)  into  (C.2.5). 


(C. 2.5) 
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Appendix  D  COMPUTER  PROGRAM  LISTING 


D.l  FORTRAN  Program  Listing  for  a  Recursive  ARMA  Spectral  Estimation 


THIS  PROGRAM  COMPUTES  AUTOREGRESSIVE  COEFFICIENTS  OF 
'HIGH  PERFORMANCE'  ARMA  MODEL  (REAL  DATA,  p=»q)  . 

DIMENSION  X(1024),EXN(30) ,EXNM1(30) ,BXN(30) ,BXNM1(30) 
$  ,FENM1(30) ,FRN(30) ,FRNM1(30) ,SN(30) ,SNM1(30) ,TN(30) 

$  ,GAM(30) ,GAM1(30) ,EYN(30) ,EYNM1(30) ,BYN(30) 

$  ,Y(1024),XA(1) ,RX(30,30) ,YX(30) ,YS(1024,30) 

$  ,WKAREA(30) ,CM(30) ,CM1(30> ,AM(30) ,AM1(30) ,BM(30) 

$  , DM(30) ,RXX(30,30) ,FEN(30) ,TNM1(30) ,BB(30) ,BM1(30) 

$  ,XS(1024,30) ,BYNM1(30) 

Nl:  TOTAL  NUMBER  OF  OBSERVATION  DATA 
IP:  ORDER  OF  DENOMINATOR  COEFFICIENT 


Nl=64 

IP-4 

NP*N1-IP 

N-Nl-1 

IP1-IP+1 

IPM1-IP-1 

GENERATE  DATA  TO  BE  MODELED 

DSEED-12345 

CALL  KAVEH(Y,N1,DSEED) 

DO  25  1*1, NP 
25  X(I)*Y(I+IP) 

WRITE  (6,101)  (Y(I) , 1=1, NP) 
WRITE  (6,101)  (X(I) , 1=1, NP) 
N1=NP 
N=N1-1 

INITIALIZATION  FOR  TIME  UPDATE 

EXNM1(1)=0. 

EYNMl(l) =0. 

BXNM1(1)=0. 

BYNM1(1)=0 . 

FENM1(1)=0.0 
FRNM1  ( 1)  *0 . 0 
A12-0.0 
A2 1-0.0 
A22-0.0 
DO  1  I-1,IP 


non  non  nnoonno  non  noon  non 
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SNM1(I)=0.0 
1  TNM1(I)=0.0 

UPDATE  PARAMETERS  FROM  IT=1  TO  IT-N1 

DO  2  IT=1,N1 

ITM1=IT-1 

AIT-IT 

WRITE (6, 102)  IT 
102  FORMAT  (/,3X,'N=',I3) 

INITIALIZATION  FOR  ORDER  UPDATE 

EXN(1)=X(IT) 

BXN(1)=X(IT) 

EYN(1)=Y(IT) 

BYN(1)=Y(IT) 

DO  20  1=1, IP1 
SN  ( I)  =0 . 0 
TN(I)=0.0 
GAM1(I)=0.0 
20  GAM(I)»0.0 

UPDATE  FEN(l)  AND  FRN(l) 

FEN ( 1) =FENM1 ( 1) +X ( IT ) * Y ( IT ) 

FRN ( 1) =FRNM1 ( 1) +X (IT) *Y ( IT) 

M=IP 

IF(ITMl.LT.IP)  M-ITM1 
Ml-Mfl 

IF(IT.EQ.l)  GO  TO  109 
ORDER  UPDATE 
DO  3  1=1, M 

UPDATE  GAM(I+1)  AND  PARTIAL  CORRELATION  COEFFICIENT 
SN(I)  AND  TN (I) 

GAM ( 1+1) =GAM( I) +BXNM1 ( I ) *BYNM1 ( I ) /FRNM1 ( I ) 

SN ( I) =SNM1 ( I) +BYNM1 (I) *EXN ( I) / ( 1 . 0-GAM( I) ) 

TN ( I) =TNM1 ( I)+EYN ( I) *BXNM1 ( I) / ( 1 . 0-GAM( I) ) 

UPDATE  FORWARD  ERRORS  EXN(I)  AND  EYN(I) 

EXN ( 1+1) =EXN ( I) - (SN ( I ) /FRNM1 ( I) ) *BXNM1 (I) 

EYN ( 1+1) -EYN ( I) - (TN ( I ) /FRNM1 ( I) ) *BYNM1  ( I ) 

UPDATE  BACKWARD  ERRORS  BXN(I)  AND  BYN(I) 


i 
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FEN  ( 1+1)  =FEN  ( I)  -  SN  ( I )  *TN  ( 1 )  /FRNM1  ( I) 
FRN(I+1)=FRNM1(I)-SN(I)*TN(1)/FEN(I) 

3  CONTINUE 

109  IF(IT.EQ.l)  Ml=l 
PRINT  OUT  AT  EACH  NEW  DATA  POINT 
WRITE (6, 100) 

100  F0RMAT(/,2X,  ’EXN(I) \3X, 'EYN(I)  '  ,3X, 'BXN(I)  ' 

$  ,3X, 'BYN(I) ' ,2X, 'FEN(I) ' ,2X, 'FSN(I) ' ,2X, 'SN(I) ' 
$  ,2X,  'TN(I) ' ,2X, 'GAM(I) ') 

DO  5  1=1, Ml 

WRITE (6, 101)  EXN(I) ,EYN(I) ,BXN(I) ,BYN(I) ,FEN(I) 

$  ,FRN(I) ,SN(I) ,TN(I) ,GAM(I) 

101  FORMAT  (2X.10F8.3) 

IF (IT.EQ.N1)  GO  TO  5 

READY  FOR  NEXT  DATA  POINT 

EXNM1(I)=EXN(I) 

EYNM1 (I) =EYN (I) 

BXNMl(I) =BXN (I) 

BYNM1(I)=BYN(I) 

FENM1(I)=FEN(I) 

FFNMl(I) =FRN (I) 

SNMl(I)=SN(r) 

TNM1(I)=TN(I) 

GAM1(I)=GAM(I) 

5  CONTINUE 

IF (IT.EQ. 1)  GO  TO  2 
A12=A12+Y(IT)*X(IT-1) 

A2 1=A2 1+Y ( IT- 1) *X ( IT ) 

A22=A22+Y (IT-1) *X(IT-1) 

2  CONTINUE 

FIND  AUTOREGRESSIVE  COEFFICIENTS  FROM  PARTIAL 
CORRELATION  COEFFICIENTS 

A11=FEN ( 1) 

DET=A11*A22-A21*A12 

CM( 1) = (A22*Y (N 1) -A12*Y (N 1-1) ) /DET 

CM(2)=(-A21*Y(N1)+A11*Y (Nl-1) ) /DET 

AM(1)=-A21/A22 

BM(1)=-A12/A11 

IF(IP.EQ.l)  GO  TO  23 

RM=X(N1)*BM(1)+X(N1-1) 

GM-X(Nl) *CM(1)+X(N1-1) *CM(2) 

ETM=1.+RM*CM(2) /(1-GM) 

DO  13  I0RD=1,IPM1 

IORDl=IORD+l 

IORD2-IORD+2 


! 
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UPDATE  AUXILIARY  VECTOR  DM(I) 

DO  14  I=1,IlRD 

14  DM(I)=(BM(I)+RM*CM(I) / (l. -GM) ) /ETM 
TEMP=SN ( IORD+1) /FRNM1( IORD+1) 

UPDATE  FORWARD  VECTOR  AMI (I) 

DO  15  I=1,I0RD 

15  AMI (I) =AM(I ) -TEMP*DM ( I ) 

AMI  (IORD+1)— TEMP 

UPDATE  BACKWARD  VECTOR  BM1(I) 

TEMP=TN (IORD+1) /FEN (IORD+1) 
BM1(1)— TEMP 
DO  16  1=2 , IORDl 

16  BM1(I)=DM(I-1)-TEMP*AM(I-1) 

UPDATE  AUXILIARY  VECTOR  CM1(I) 

TEMP=BYN (I0RD2) /FRN (I0RD2) 

DO  17  1=1, IORDl 

17  CM1(I)=CM(I)+TEMP*BM1(I) 
CM1(I0RD1+1)*TEMP 
SUM=X(N1-I0RD-1) 

SUM1=X(N1)*CM1(1) 

DO  18  1=1, IORDl 

SUM= SUM+X (N 1+1-1 ) *BM1 ( I ) 

18  SUM1=SUM1+X (Nl-I) *CM1 (1+1) 

RM1=SUM 

GM1=SUM1 

ETM1=1.+(RM1( 1. -GM1) ) *CM1 (I0RD2) 

SET  VECTORS  FOR  NEXT  ITERATION 

DO  19  1=1, IORDl 
AM(I)=AM1(I) 

BM(I)=BM1(I) 

19  CM(I)=CM1(I) 

CM(IORD2)=CMl(IORD2) 

RM-RM1 

GM-GM1 
ETM-ETM1 
13  CONTINUE 
23  CONTINUE 

PRINT  OUT  AUTOREGRESSIVE  COEFFICIENT 
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WRITE  (6,105)  (AM(I) , 1*1, IP) 

105  FORMAT (/,3X, '  RECURSIVE  SOLUTION  -  ' .//.10F10.5) 
RETURN 
END 


NOTE:  Above  program  may  be  applicable  to  complex  data  by  making 
following  changes 

(i)  Declare  all  variables  to  be  complex  value  except 
integer  variables  (i.e.  IMPLICIT  Statement) 

(ii)  In  DO  loop  25,  take  complex  conjugate  on  the 
variable  Y(I+IP)  (i.e.  Y(I+IP)«CONJG(Y(I+IP))) 


nnnn  non  oonoon  ooono 
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D.2  FORTRAN  Program  Listing  for  Generalized  Levinson's  Approach 
of  ARMA  Model 


Generalized  Levinson's  approach  discussed  in  Section  5.3  is 
programmed  for  the  premodified  method. 


THIS  PROGRAM  COMPUTES  DENOMINATOR  COEFFICIENTS  OF 
'  HIGH  PERFORMANCE  ’  ARMA  SPECTRAL  ESTIMATION 
BY  GENERALIZED  LEVINSON'S  APPROACH  (REAL  DATA,  p-q  ) . 

DIMENSION  X(64) ,FEN(30) ,YV(30) ,XV(30) 

,FRN(30) ,FRNM1(30) ,SN(30) ,TN(30) 

,  BMN  (10 , 10)  , BBMN  ( 10 , 10) 

,Y(64) ,RX(30,30) 

,WKAREA(30) ,CM(30) ,CM1(30) ,AM(30) ,AM1(30) 

,DM(30) ,RXX(30,30) ,BM(30) ,BM1(30) 

Nl:  NUMBER  OF  TOTAL  OBSERVATION 

Nl=64 

IP:  ORDER  OF  DENOMINATOR  COEFFICIENTS 
IP=4 

NP-N1-IP 

N*N1-1 

IP1-IP+1 

IPMl-IP-1 

GENERATE  DATA  TO  BE  MODELED 

DSEED- 12345 

CALL  KAVEH(Y,N1,DSEED) 

DO  25  I-l.NP 
25  X(I)-Y(I+IP) 

WRITE (6, 101)  (Y(I),I-1,NP) 

WRITE (6, 101)  (X(I) , 1=1, NP) 

101  FORMAT (2X.10F8, 3) 

Nl-NP 

N-Nl-1 

INITIALIZATION  BASED  ON  THE  FIRST  TWO  DATA  SAMPLES 
X(l)  AND  Y(l) 

DO  40  1*1, IP1 
DO  40  J-l.IPl 
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BBMN(I,J)=0.0 
BMN(I,J)=0.0 
40  RXX(I,J)«0.0 

A11-Y(1)*X(1)+Y(2)*X(2) 

A12«Y(2)*X(1) 

A21=Y(1) *X(2) 

A22-Y(1)*X(1) 

RXX(1, 1)=A11 
RXX(1,2)=A12 
RXX(2, 1)=A21 
RXX(2,2)«=A22 
BMN(1,1)=-A12/A11 

FRNM1(2)=RXX(2,2)+RXX(2,1)*BMN(1,1) 

C 

C  SOLVE  FOR  DENOMINATOR  COEFFICIENTS  (AM(I) ,1=1, IP) 

C  AT  EACH  NEW  DATA  POINT  FROM  IT=3  TO  IT-N1 
C 

DO  38  IT=3,N1 
IPM1=IP-1 

IF(IT.LE.IP)  IPMl=IT-2 
C 

C  UPDATE  ROW  VECTORS 

C 

DO  37  I-l.IPl 

YV(I)*=0.0 

XV(I)=0.0 

IF(I.LE.IT)  YV(I)*Y(IT+1-I) 

IF(I.LE.IT)  XV(I)=X(IT+1-I) 

37  CONTINUE 

DO  39  I-l.IPl 
JF=I 

IF(I.EQ.l)  JF=*IP1 
DO  39  J=l, JF 

39  RXX(I,J)=eRXX(I,J)+YV(I)*XV(J) 

A11-RXXC1,1) 

A12*RXX(1,2) 

A21*RXX(2,1) 

A22=RXX(2,2) 

AM(1)=*- A21/A22 
BM(1)— A12/A11 
IF(IP.EQ.l)  GO  TO  23 
DO  13  IORD-l.IPMl 
IORDl=IORD+l 
IORD2-IORIH-2 
C 

C  COMPUTE  AUXILIARY  PARAMETERS  FEN (IORD+1) 

C  AND  FRN(IOREH-l) 

C 

SUM“RXX(1, 1) 

DO  27  I-l.IORD 
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27  SUM-SUMfRXX(l,I+l)*AM(I) 

FEN  (IORD+1) “SUM 
SUM-RXX(I0RD1,I0RD1) 

DO  28  1=1, IORD 

28  SUM=SUM+RXX(IORDl,I)*BM(I) 

FRN (IORD+l)=SUM 

C 

C  COMPUTE  PARTIAL  CORRELATION  SN(I) 

C 

SUM«RXX(IORD2,l) 

DO  29  1=1, IORD 

29  SUM=SUM+RXX(IORD2,I+l)*AM(I) 

SN (IORD+1) =SUM 

DO  14  1=1, IORD 

14  DM(I)=BMN (I , IORD) 

C 

C  COMPUTE  PARTIAL  CORRELATION  TN(I) 

C 

SUM=RXX ( 1 , IORD2 ) 

DO  30  1=1, IORD 

30  SUM=SUMERXX(1,I+1)*DM(I) 

TN (IORD+1) =SUM 

C 

C  UPDATE  VECTOR  AMI (I)  ;  FORWARD  SOLUTION 
C 

TEMP=SN (IORD+1)  /FRNMl(IORW-l) 

DO  15  1=1, IORD 

15  AMI ( I ) =AM( I ) -TEMP*DM ( I ) 

AMI  (IORDfl) —TEMP 

C 

C  UPDATE  VECTOR  BM1(I)  ;  BACKWARD  SOLUTION 
C 

TEMP=TN  (IORIH-1)  /FEN  (IORD4-1) 

BM1(1)=-TEMP 

DO  16  1=2, IORD 1 

16  BM1(I)=DM(I-1)-TEMP*AM(I-1) 
SUM=RXX(IORD2 , I0RD2) 

C 

C  COMPUTE  AUXILIARY  PARAMETER  FRN (I0RD2) 

C 

DO  31  I=1,I0RD1 

31  SUM=SUM+RXX(IORD2 ,I)*BM1(I) 

FRN (I0RL2) =SUM 

C 

C  SET  FOR  NEXT  DATA  POINT 
C 

DO  19  1=1, IORD1 
AM(I) “AMI (I) 

BBMN(I,I0RD1)=BM1(I) 

19  BM(I)=BM1(I) 


u  o  o 


13  CONTINUE 

BMN (1,1) =-A12 /All 
DO  43  I=l,IORDl 
DO  43  J=2 ,IORDl 
43  BMN(I,J)=BBMN(I,J) 

DO  41  I=l,IORD2 
41  FRNM1(I)=FRN(I) 

38  CONTINUE 
23  CONTINUE 

PRINT  OUT  RESULTED  DENOMINATOR  COEFFICIENTS 

WRITE (6, 105)  (AM(I) ,1*1, IP) 

105  FORMAT ( / , 3X , '  GENERALISED  LEVINSON  SOLUTION  «  ' 
,// , 10F10.5) 

STOP 

END 


NOTE :  Above  program  may  be  applicable  to  complex  data  by  making 

changes  as  described  in  Section  D.l  (See  Expression  (5. 3. Id)). 
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